{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Nonlinear Waveguide Demo\n",
    "\n",
    "This is a demo of `angler`'s nonlinear solver. \n",
    "\n",
    "In this notebook, we'll make a straight waveguide and give it some Kerr nonlinearity.\n",
    "\n",
    "Then, we'll simulate the waveguide and measure the self-phase shift as a function of input power"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pylab as plt\n",
    "import numpy as np\n",
    "import scipy.sparse as sp\n",
    "\n",
    "from angler import Simulation\n",
    "from angler.constants import *\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "omega = 2*np.pi*200e12         # angular frequency (rad/sec)\n",
    "dl = 0.01                      # grid size (microns)\n",
    "Nx, Ny = 600, 200              # x and y grid size\n",
    "eps_r = np.ones((Nx, Ny))    # relative permittivity\n",
    "eps_m = 12.25                  # maximum relative permittivity (Si)\n",
    "wg_width = 20                  # waveguide width in pixels\n",
    "eps_r[:,Ny//2:Ny//2+wg_width] = eps_m        # define a waveguide\n",
    "\n",
    "# define the region where nonlinearity will occur\n",
    "nl_region = np.zeros(eps_r.shape)\n",
    "nl_region[100:Nx-100, Ny//2:Ny//2+wg_width] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "linear power transmission of 0.9996562161761373\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEoCAYAAAAqrOTwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvVmsZEl6mPdFxNlyvXstXdU9izkcLqJGHDZnRFEcU7Apz4s8LzQ4JmBTgATDD4SfBFN6kSX6RbJfLMAEbEKi4QU2DVCQ0AJo0DJkitSI5EwPyTHJ4ZDT3ZzpruqqunXXvLmcPEuEHyLOyXPz7ln3VmXVxAcU6t7MvPmf/4+I/4/4YxPGGDwej8fjedWQL/oBPB6Px+O5CXyA83g8Hs8riQ9wHo/H43kl8QHO4/F4PK8kPsB5PB6P55XEBziPx+PxvJL4AOfxeDyeVxIf4Dwej8fzSuIDnMfj8XheSYIX/QAej8fjORvZv28o0mv9TjPZ/TVjzOev9UuXEB/gPB6PZ5kpUoJP/ofX+pX57/9Pm+e9L4T4PPCPAAX8Y2PMP5h7/3PAfwf8eeCLxphfabxXAn/gfn3fGHO9D38FfIDzeDyeZUYIhFTPUZxQwC8APwE8AL4ihHjLGPP1xsfeB/468LdO+YqJMeYv3PiDXgIf4Dwej2fJeZ4BDvgM8I4x5j0AIcQvA18A6gBnjPmWe08/zwe7Kj7AeTwez1JzIyO4TSHE243ff9EY84vu53vAB433HgCfvcJ3J+67C+AfGGP++bM96uL4AOfxeDzLzM2kKHeMMW9e95c6PmKMeSiE+Djwr4QQf2CMefeGZJ2LD3Aej8ezxAhAqOeaonwIvN74/b577VIYYx66/98TQvw68IOAD3Aej8fjmUMI5POdg/sK8AkhxMewge2LwE9f5g+FEGvA2BgzFUJsAj8K/Dc39qQX4AOcx+PxLDnPc5GJMaYQQvws8GvYbQK/ZIz5IyHEzwNvG2PeEkL8MPDPgDXgrwkh/r4x5vuB7wX+R7f4RGLn4L5+hqgbxwc4j8fjWWae8zYBAGPMrwK/Ovfa3238/BVs6nL+7/4t8AM3/oCXxAc4j8fjWWIEIKQ/VXERvNU8Ho/H80riR3Aej8ez1Dz/FOWrgg9wHo/Hs8y8gDm4VwUf4Dwej2fJ8QFuMXyA83g8nmVGiOe90fuVwQc4j8fjWWLsKkof4BbBBziPx+NZZvwc3ML4AOfxeDxLzXM/quuVwQc4j8fjWWaET1Euig9wHo/Hs8QIvw9uYXyA83g8niXHB7jF8AHO4/F4lhm/yGRhfIDzeDyepcYHuEXxAc7j8XiWGfHcb/R+ZfC3CXg8Ho/nlcSP4Dwej2eJ8asoF8cHOI/H41lm/CKThfEBzuPxeJYcH+AWwwc4j8fjWXKkFC/6EV5KfIDzeDyeJUYIgfABbiF8gPN4PJ4lRwgf4BbBBziPx+NZcnyKcjF8gPN4PJ5lRuBTlAviA5zH4/EsMfZGbx/gFsEHOI/H41lqBNLPwS2ED3Aej8ezzPgU5cL4AOfxeDxLjg9wi+EDnMfj8SwxQvhVlIviA5zH4/EsOcLf+7IQ3mwej8fjeSXxIziPx+NZcvxJJovhA5zH4/EsMUIIPwe3ID7AeTwez5LjV1Euhg9wHo/Hs+T4ALcYfpGJx+PxLDMCpBDX+u9CkUJ8XgjxJ0KId4QQf/uU9z8nhPhdIUQhhPjJufd+RgjxTffvZ67RElfGj+A8Ho9niXneZ1EKIRTwC8BPAA+Arwgh3jLGfL3xsfeBvw78rbm/XQf+K+BNwABfdX+7/zyefR4f4Dwej2epee4Xnn4GeMcY8x6AEOKXgS8AdYAzxnzLvafn/vY/AP6lMWbPvf8vgc8D/8fNP/ZJfIDzeDyeZeb5n2RyD/ig8fsD4LPP8Lf3rum5rowPcB6Px7Pk3MA+uE0hxNuN33/RGPOL1y3kReMDnMfj8Swxdg7u2r92xxjz5hnvPQReb/x+3712GR4CPz73t79+1Ye7LvwqSo/H41lmXIryOv9dwFeATwghPiaEiIAvAm9d8ml/DfirQog1IcQa8Ffday8EP4LzeDyeJed5LjIxxhRCiJ/FBiYF/JIx5o+EED8PvG2MeUsI8cPAPwPWgL8mhPj7xpjvN8bsCSH+a2yQBPj5asHJi8AHOI/H41lqxHM/i9IY86vAr8699ncbP38Fm3487W9/CfilG33AS+JTlB6Px+N5JfEjOI/H41li/IWni+MDnMfj8Sw5/izKxfABzuPxeJYYIUD5ALcQPsB5PB7PkuMD3GL4AOfxeDxLjED4ALcgPsB5PB7PMuNTlAvjA5zH4/EsMQIf4BbFBziPx+NZYoSAwAe4hfABzuPxeJYYP4JbHB/gPB6PZ5kRfpHJovgA5/F4PEuMHcH5UxUXwQc4j8fjWXL8CG4xfIDzeDyeJcafZLI4PsB5PB7PEuM3ei+OT+x6PB6P55XEj+A8Ho9nyVHP+cLTVwUf4Dwej2eJ8XNwi+MDnMfj8Sw5PsAthg9wHo/Hs8T4o7oWxwc4j8fjWWL8KsrF8QHO4/F4lhwf4BbDBziPx+NZYvwik8XxAc7j8XiWGH+bwOL4AOfxeDzLjB/BLYwPcB6Px7PE+EUmi+MDnMfj8Sw5PsAthg9wHo/Hs8T4RSaL4wOcx+PxLDF+kcniXCnAbW5umo+88Qbmpp7mOXJT1eUqtrnsM9yUvZ/VBhc91zLYuMlZz3PT9fllck3XaYum3jdh4+u262nPuIiMb7//Pjs7Oy9Tsb+yXCnAvfHGG3zpS19a2gD3stSoZbXfWdy0Xc+yx00H4GXjedt5We17medaRPZN2Lf5HNX3/+iP/uj1CvEpyoVZKEUpWA7n8TIU+TLYaRGep22vQ5a3883Lel42vi45z8O2z0eG8NflLMhLOwd3U431OqvRy+h0X7Zm9DLaGF4uO79sNr7pkdpNyjkL6QPcQryUAe5l6Yl6Xj5uet7IY3mZ7Pyin08Ayse3hXgpA9xNYni5etee6+G0Mq9eu24Hd13163lkIa6b52nnZ+Wi53luvkKA9HNwC/EdF+CWrRFdhdOq+DLr8zI64NNYljnnJuc9T/XeMtn5Ms+yjHa+iOcR5OwIbplK8+VBvugHuCrPo5iXsZGdpbc4570XheFyDnhZuKzzXRY7X9Z+y2bny/Cy2fiqn10UKcS1/rsIIcTnhRB/IoR4Rwjxt095PxZC/J/u/d8RQnzUvf5RIcRECPH77t//cO3GuALXNoJ7UXuMboplSlUuy3NcJ8ti32V4hlcdb+Nn43nPwQkhFPALwE8AD4CvCCHeMsZ8vfGxvwHsG2O+SwjxReAfAj/l3nvXGPMXnt8Tn821jOBelgr8Mgbby9p2WcrgZbTxq4wvj6uzdDYTAimv998FfAZ4xxjznjEmA34Z+MLcZ74A/M/u518B/j0hli+P+swB7iKNlk5jz9KwdI7kEvj6fPN4Gx9HcCMpyk0hxNuNf/9ZQ+Q94IPG7w/ca5z2GWNMARwCG+69jwkhfk8I8a+FED92Eza5LM+UorzK6OI6nJmv+J7vdF7GTsF3Ajedcr+BFOWOMebNa/9WeAS8YYzZFUL8EPDPhRDfb4wZ3ICsC1k4wL1swWYRx/Ci54muKvtFr0J7GW3s8Sw71QjuOfIQeL3x+3332mmfeSCECIAVYNcYY4ApgDHmq0KId4HvBt6+8ac+hee2TeBFO1+PZ57vlMD6IjsRi8r1/qLB8z+L8ivAJ4QQH8MGsi8CPz33mbeAnwF+C/hJ4F8ZY4wQYgvYM8aUQoiPA58A3nt+j36c77h9cB6Px/My8bxHcMaYQgjxs8CvAQr4JWPMHwkhfh542xjzFvBPgP9VCPEOsIcNggCfA35eCJEDGvjPjTF7z+3h53iuAc73yjwej2f5Mcb8KvCrc6/93cbPKfAfnfJ3/xT4pzf+gJfEj+A8niviO2qe540/i3IxfIDzeDyeJUZwudNHPCe5coDzq948Hs+rylKOzP2FpwvjR3CvGD595vG8WthFJi/6KV5OfIDzeDyeJcffJrAY1xvgjJ79LK52CphuDDuu1FtpynwGuVfuIT2DXD03xFpY32eQeSW512TjK8l8RrnfqTZ+JvteUe4y2PilqE/PyAvY6P3KcH0Bbr7yVL9fUIlOcwqXbqynNdAF5T6TzOr1BXStXn+eur5IuQvLvKTcZdL1UnKvWddLdyZesnp8nky4OV3Pkn0pXa8LAeqlu9hsObieAHdWBareO6MSnVVpm++fWYnOk/kMcp+p4l6jrsce4YZ0PU3uie9dQO5VZR5buHSRzAu+9ypyj/GibHwRl3TCV/7OBWW+KBtfxMI2viC43pivuAJ+BLc4L+cc3DM4wRuX+yIa6UvIi3L6123jixz+udxgMK8+s1SBdUGW2cbPB+Hn4Bbk2QNcVYHOy+ef0lialee0eiQan7tyj3BBuefKbMq6SO4cF+nalH2Ci+TelI2fQe6xrzlH7skP35yNz3UPL0LX855hXu55f/YCZF4k90XVp2fyFeewsMxrwo/gFmehAHdiL9yCOe7zOkgX7rdbMC16ntzqdcE5FfcG5FbvLaTvM8rkNLmXCTQXpCcvo+sz2fgM+Zcp23O/+zyZc1ym43Kmrud1DOflzsm/yMbn1uOL5DZlLSD3Qsd/g/XpTFnP28bXjZ+DW5hnG8Fdpld2wcjitN/BVphTK+5le/kXyD0r/XCpyeoFRo6VuIvkntv7Pa+hzul8VV3PlHtFG8PldF3YKS0o90xd5+VdVI/nZDdFXVSPzww2l5E7x2V0XcjGlczq5zk5F8k9t81e1sbzf974+UpyL5J5SbkL+YprxI/gFud6UpSu8oi5SnSiXsw5hrrRmOOfFELUzqD63AnHcJbMZgN9BrknRqhzDeVcuY1Go83JhnKjcpvvnSMTQDtpJ+TOyT5T5jmyL6Nr9bmzAvq5cs9AmwV1df/PywRXfhfIPk/uZTpN58o9g7Ns3NT1PLlnyryE3PN0vSjIXWjjU+Rftj6dyhXr02mdw0u32RvCx7fFuLaB72mV9rTXLoMxBmPMhZO8x75fn+0o4GQv7LQG2nytGRzOlKv1hXLn5Vcymp8+7VmOP5h1CovauPn9Vy2Rq9gYLg6qZ+paBdR5Z6Rnss4r2/PK61y5Dfnn2vKc9+Yd/vwnm52c07hQ7nn2voqNL/F9ly3bivl6fG6bvUjuBXX5NBvX7emCOnBu25mre+fJPe81z/JxLasohXGOvllBXK9IyNN7wE2noIFmfal6K9KY07sulTPSjQoKUNrvrzOFZ+T0jTGNBnJS5pm97qbcY45/pvJZPeDTdC0vKfeYrk19T7Mxs0Z+lq6VXFk/k5idVt5o7BfZeF7XY3NSl7TxaWY+IbfUJ3VtyG1yllxZvydOTROeqquTIeTpvf3KoV5kY3FGHT7Txk1d59LeTUd+5Xrc0OtSbbbx2rFHd3LPqscIcTygz3dcmnKr75en97XP01U3ZF6Ljc+Re5Z/unCEfk1IfwLwQjz7HFyzwjYbS+UMNHaGtOEc6pRd9VEDulGDqsLUAsQp6UlxlsyK6iVx0iE1K2xTJsbmubWYVVx1WbkNXU9rMNqcreu83BMBfd4hnKXvOTY+IRNr48pBqMbnj8m4wManOf55p3BVXYXRUBYndb2EjY0xlKfoigEjxCyon2fjeV2NBiMRyjYV0wg09vczdG3IVcw5wtOc/Xm6won5sKvYWBuOdWCuZGP3WtPhn2XjSm7VZpty69FSWZyUCbWN53U918ZOZgm1jdU12LiScBkbK24egU9RLsozj+BsBSoQxpx0DsJgpP3MvEOyaUhbceZTCwZjJ1YRGOYcUqPCCl3OXmvKdIFNGF07JJgFGmOgPCWdYjAoZs7BVHKbwe0CXdHYwAq1YzhP13m54jRnVPdEzQldgQttfJauwjXUEzZmljISujxdV6GB4PRgc4GNhQGkON55aTreU3WVYMpaLsI+n2Zm4/IMG1tna+qgLox1Ggp9CRtLDO65hJyVLcfr03lyhRAnA2ujPlmTnizbWtdTltBpoNSn1ycbWWc2PnGXmL5A10pAo4PYDKoX2dicFmz05WwsVGDbLKCRtY1LJ+wsXwG2I2HOsLGoAtS8jV2bPTZqhWNyz2uzla6Cm72z7Ttlf+x1c6UAd+rEsXYVVs/1zABkYP1Q5RgajaVsOMD5HlLV2zZY51A30sYoRuiydg5170xIl1IKbIPRunZIs16ZlVvqk71QY7v4dbA51jtrOoVK12bKSthxpzkl09J0ChfJpdnznQ/m873QKs1S6TfXy9fMbNx0EGLOxnDcxlWgqW18IsBJhAvedbBx322MqZ3CWboqacvhhCM8T9dKrrQOUBh93CE5G5cX1Cc112kS8zau6tS8rs1yNRqQxzouZcPG83Kr1QjNYHOsAwFnywWElHVHrWnj03QFjq3gkZJZYK07asWZNhYyON4xpZHyhjPlllRXuswF9Gbm4ZyyFTJwwUYfWxlQ2bh6jjNt7OQ222xt47P8U9Vmq46EnI3Om3LPa7Ondl6uG+FHcItyfaso68rUGEnoAmRwvJKfSLUc7+k3e4II+34zsoo5eXWjpdEDMxqBrHuC9cvGoJ2s+V5o3ds2AoFBudQOVaA5S25TX6PBiGM2qHSsZF0k1+p7vBda93jn9XUB3dp31rSbwbxyCid7oW5Ewxk9/LmyPSazRp1wGprL63ra6Qzn6tqw8WkYZvWJhtx6hAyny52zcXNBghHSdtSMHWWcCKxOx2YZz57HBdQqvhiDNgI1N6I5z8bCSMyc06+e7rSRoxS2vKVwnSU49rx1h7QKdPO66gJUdKJcqznkpo3nR3DG2KI5q82eJrdZrsK4NnvCT5gT7aiSWY8azwoAF9m48hWn/KludITPqsfzut4EAuHn4BbkmQJcHWzKDKGLOq8vwFZQl1MXzgmbRu+onOvpV04J7NUQQthKLZRrMFUBN0aMosxozk8JsDLr0Vx8rBdaNfZSGwoNugoo4BYiuN69m0BQai7Q6HKma1OukOB62YLIOqRqVOOCaqXjvK6VXDBIJY6tSjthXzeyqt+v5oWaNm6MLiobV6PlY9MIQqBklaYE1bRxYzRTl+u8jWVgZQrrtCobGzPr8c7buCrbysbS2LkT25Oxjq8u12YQOMvGjR53qa1tc21q+zbLVgs3sJGNXn5TRiV3bo5ICAlKu5FUYOUqV48M9cjtNBsLm2cHbUepuumjmjauynVeV+map1Sz50W6ANPIRMzV48rGwlgb139r9HE7N9rsTNdgpmvpOqhCzkZvDRvP62r1Ewg3sqnfnm+z87pCbeP5oYp2di60y0SYk21WCFuH6s5LI/txYZsF5ysEmGqe1czK9qI263S+6TMP/QhuMRYvl0YaSxQZosyPp1mEBF1C4EYzOnRpysopUDvAorQVqkIJ20AiJSjNXKqw4YxEkc8abFOm0qAC0M4hQZ1yyEvreOcbqXUMbuoBQAr7Nwhk0xkV2cnUnevlE7ieoQlngRWrbyWvKI+vyqrk4hYhHBtNuQYqdAGV3GZ6R5d1T9cICaUCFTXmDzhV18rGgesZSgFautRoU9fKxvPpHWdjA1ZvHdZvlXVgPd3GJVZuNWdRrX4TzVRSkZ1MK1WjVVXVJwVEx+bfcm0otA3qTRsrAUa6jhMGI5vO19XhMrNyy/y4rrLRYYJZsKHq4ZvaAed63sb2GSKFq8ezDkTdcTmrPlW6BiD07BkqJ1/JO6s+VdN2SoA2ou671GWbpydTslX7Cd2oT7oRejU6Yya70rv+UwFhla0XtoMYSnFc12Y9brZZUXVcqDMSxnWISyersm3TxpWuVUZWCOFGmlVnuJi12apc59ss2M6ayzRVbbYq2/NsXLVZIWczgTeFn4NbjIUCXG3rOtjkiCKlmsytJ3pVZD8mA/uerNJ2tuJUDjArZ+kOgW2clQwpsA3FyauckcinbgRX1PMYNi/uUqK6tI5QWedrA41tlFmpyVxvtGqjVq4hMqLunekq2lQT1WV2QtfjcpPaSVQpEGOqEaNtMKfJDaR10kJIOz/VdISVrkXaSO84uVLVjRQhMU7X2RyNtXFWWvkwkysFRMycYalx80SNHn7hdC2zU3R1HYeA2nEYI+rAVpSGXFvHf9zGtiwkdi6uHmHoStcUUUxdj3tO1+aS8jrwzOaGCqdr5QSbuobGBnWQrjcuZmmzyvEWtk7V82KVXFePCaJarka6wEZt49PSdqGZjaqajrDquNS6zttYSFufcA7YOV/T6DxU5XpafVLuhcB11BDM5t+KWT0+oWtVnwJtdWUW2Jpyc316ihJXp5Bi1qabnbQ5G9e6qqjRYbIpWe3KsNDUneCm3KrNVilnJZi12aavcDKbc55GiFkqVmnno6r2Q12259kYV5/FmfnR68PHt8V4thRlHWzGmMkQk6Wz96RCJB27xFUqCBKM1hg1SzlMS02pYVKUTAvrMKUQhEoQSglYhx9VcxdgK2SRWcebjjBFji6yWiZBiGh17bxJJZdZ6iwrbZCZ5Jq00PXksRSCJJDoQNYpj7Ja7dEYzYh8DNMJJkvRrrEIqRBRgoi1CzRRnVOoUnWVU5gWp8tthRKlDWUjj1WPjosUpmNMkUGRz+RGCSKI3JyfQASJc7428Vc5o7TQTAtN3thLFwcSYxRKGqSQBHIumOep1TVLMemollnLTdqzgONsXDmgyimcZuNQCdqhQqIJpKxTPNb5WIcv8wkmz47bOIgQQdiwsU3tVenJysbT0jDOShdsrNxQStqRsr1tYYj0zPlSZrUTlPkEPRnV9amSK1sdmzo2xpYts05EVmpyzZk2bgXKBThzbLRad1zyMRSF1XW+HidVJyJyTt+WbdUpnBaGSVGSl8d1jQNbnhIopOuoVSm7yuHnY0w6RjfbbBDZsnXtzKionrNqBrfKxk1dK7naSBdwjBtNVfVpCmV2wsZCKpBqZmMVYHRse07CZlxybUhdJ7gpt2o7oRuu2mDOzMZOrihSzGQIWh+Ta+txxwW5ABPE9Qi5KtvTbCyFrcMGgZK2O1zP198Awsn0XJ1n3wdXFuCckZmmNsUhFShbcUUQIHSC0UU9mqpGNTbgaMa5Ji+tI5RCkBiJlhAqQdl0RtBIKeXoLMUUOSbP7KhJKkRsnT4hx1JNupoLcwGnckiVI6wqkBSCSBmkmY2khDFWTydbZyl6MrK6Ynu+EpBhZPVsTtxDPZLKS8M4L2unL8XsGoxQ2XRMNaF+bHGLtg7QFNlxG2sNQe5sXB5bVKPrNLANquO8PKGrEoKIal7U9VCrQONGNiZL7b/cOX2Xoqs6D0IcX8xTdSTKU5x+VbahlITSpZ+bc0SlnRfS0wkU+XG5scYUGVIqCE6uiNPYsk2LkrSc1SeAOABZQBKoOoVqy1XP5JYZ2nVczHQy01VrtJQIFdpefp3mUsfmaSobN3UtjbVxIJW1b2O0Wo+UXXAz6WimaxhBkVsnrCKELmobN+1rHbB2QcC+H7ohsVZWbtxsO9VCC11ANsWk47rt2D/OrV3CCIS0crV26WtTpydPs3ElN1SC0phZOraSq2dl27SxkQoRRjMbl7NtOMat8S2NqXVMS820KGsb494PpGp0asQxuZV/ovIVlY2xbRbAlHFdp3RdrrP20+yohbIKqIEd2Z2ycvq68fFtMZ5tDs71uPXwAH24O3P6UiHCENnuI3WJUBHoVp13z12qbjAtGOclw6xknJdoY+qeYOxyAEqI2dLgyvnmY8qjffTRwbHAKoIQ0epYBxIlCBVhSluhmymHwzRnb5IzyWcBLpSScViykoQEMrQjzcqHVmnYfEx5uIs+2sdM7OgRQAQhJu1gihy5KhBRq071VKPG4dQ6hcPUys21rgNcN1I28LRCInXcvjKfYEYD9OFu7RhMkde9XoIQBYh+gChbYDSllnW6bjAtOJwWDKeFS++YugfajTTtUKFkQBI0UjtupGqGh5SHu9b5Tl1Pv5Lby1FhZB1v2XI2Nq7DYsvzMLXlOy1nusaBJC8NSkZEyrg5IlC6RBRTSEdW1yLHjAa1rrbjEoLWyCBABLHNCIgq7WwYZiX7k5yjrHQpUis3Dkp6kSIONButkETZERhauxHyFDMZ2jo8PjrWiRBxgmz3UFIhYo0I7cKlkioLYUcze2l+wsZWV1sPQmkXNmhDnYWQ+YRysIuZpujxAJPnth67TprJMys3SBrO141oCs1hass2L3Vt41AKunFAoqTrrLn2I13KLp/C1NZjMz46vc0CstWhDGZOvzQ2TTjJZ3KnjY5aJdfqapfua9zcapnVo6h5G4sghDBCtntIsG1WF6BDtMJlPYyTV3KUlbXcqs0mwWwEV1YezZWtzCeUwwPbZrN0ZuMgRHT6GF0ikw4iSOoUcZXunpaG/Ule61r5ilhJ2qFyMkOXnr3BIRzHFtF6rsCzpSjzFJlPyPa3KXcfU47HlFmBigJkGCB7I5QuCaIEEXXqBRCFts5oZ5wxyTV7k5xRVpAVmiiQ9JOA2KUd4kDScbPXwmhEPnNG5f5TTDoiH00wpUZFAarbRWUpst0naHUwuuecvnUMw2nJzjhnZ5xxNLUyAaJA0osDcjc53o+D2copF1TN8JByf9vKHh5RpDZ4BkmE6vZmgTXq1ivgctdQ9tKccV6yO84YpFaukmJOrqQVyNmcSZnBdIw+3KXYfYTJUsrhsLaxarcRUQK6RAUhRB2AOj05zq3D3xlnHDpdKxt3ooD1Vsg00iSBpB1a52ud4BjSEeX+NuX+U/RoQJFmx21cZDbYdVcRSa8OrJk2HKYFR1nB7jjnIM1rGysp6CcBeakJlSBW0WzJvi4QRUrZKNdyeESZ205R2GnZYF7k1jmFbTAaTdXLts5oe5RxmBZkpWaSlSjneCdJQCtUhFLQCgSlkS6Y29RkcbhLufsYPTyo65NQkrDTwvQ3rK6dPiLq2lGdMHVn6Sgr2R5OOZwWTLKy1rUTBUwLjRKiTkEbY1zqeep03bYpu6ODY/VJhBGqv45QChklCL2CMZrSSCt3WrA9ytib5ORa13JbkWIlL20aWNoU3rEFNUWKHuyi97cpB3uU4zHa9eTgVP3jAAAgAElEQVSCJEL2RgSAKTJE2HKZAZfqLg17zsb7ExvQS21QUtCKFOvue9qRTcuW2o7ERD61nbTRwNanwW5t4yCJkFGI6a5auyVtcG3WuKA6mBbOV9hOUyW3FSl6cWB1FYJWoCjDRpstppRH+64eb2OmaW1jFQao/tiO6topKm5B2QOOd4R3xhn7ac4gLeqOdhRIVpOQ0hiSQCKFuvIZr1dBNBZjea7GM43g6gr09CFH7z8hOxqhswIZBagwINkY0tIa2eogOmuYMKkby2FqG+jBJOfRwYSDcV47341uRDcJ0QbaoaIXSQzK9cqsEyyePmTy4Tb5KCU7GmFKjYwC4tUe3Xs5ejJCtDqIZAWwlXZaavbSnEdHUx4eTHg6SG3PTBvakWK9G3N/rUWiJNoYupELrEVWj2bK3ceMH26THhxRVo0liUhWe7SLHIKQoLuKcOnYQtuGsju2jvfB/oTtQVoHuDiQbPUT2zMMJCtJUK9OE0WOHh5Q7D4i/fBD8lHK9OCotnHU6xB2Elq6REQJsruOMDZlVTRs/MHBhMcHE6aNALfaDrm/3qYXB8SBpBe7oWOZ2eDmdJ08ekK6e0g+siO4ysadO6lL4ZWItnVOueu4bI8yDqc57+9N2BtOOUoLlBQoKdjoRtxfbxMq6UbLLl1YZtYJ7m8zffSAYpQy2T2kTLM60ITthLYLcEF3BaHLuuNinWDOw8OURwcTJlnJ2AW41XbI3dUW/SQkUdbGug40qa3Du49JP3zI5OnxABf1OrRuTYgBihzVXbXBHJiWhp1xzmGa8+19W5+O0pnTr2yshKAbK4Tbi1eNzPPDXYrth2RHYybb+3UnIuwkqCS2Nla2E0EnQ5gWpRu97YwzPjxKebA3ZpKVHKW2I9BLArb6CavtECkFK3FAXkqEcKnu8YByf5vs8QPG2wd1m61snGyMSQCVpYS9NbQu0FiHbztoOR8OrNyDcV7r2ksC7q7a5+tFCpGEth6bWZstD3eZPnrAxMk1pUYlESoKaW2NSMB2IlprmLBlU85uFPX4aMqh8xWVjduRYquf0EusG+tFASuxmi0wcRmX8ulDhg+fUqbZrDOcRLQ2hrRey1DTCbK3iuhsYHCpX1e2j4ZTHuyN2R1m9YrOysaladGLAiI3vXCTR5n4VZSLsXCAE1XF3X/K4M8ecfjuQ8a7I4pJSdBSqFDRubuKKTXdpE24cRdRttBYx7A9ynj36Yi94ZT3no4YH03RpSYIFZ1uzEo7pNSGbqS41XHL0HVh0yu7jxl++0MGf/aI7GhMup9SZpqgpejc6lGmGclGH9VbI9h4DbCVdpJrtodT3tk+4puPhxwMUspCo0tNGAf0ujGTrCAO7Eq7212bnxfFtB7NDN59wNEHTxg9OSIf2RRl2Anp3O5R5gX9KEFt3EEECdqlznbGOe+7wPbNx0MGg5QiL5FKogLJ9mqLrCiJleR2J0Ib10jzMcX+NumHH3Lwpx+QDca1jVUkaW+2CdoJq1lBN+kQ3LpPWRaUJiQtDTvjjG/tjfnGhwO29ycUeVnLbfdihmlhHUQUcKcb29FFkdk07P424wePOHj3IaNHB2SjHFMagpaivdGhGE1YDSPQJeHGXTCaQkv2JzkfHqVsD6Z848MBe4cp+dQ6XxVIOt2YSWbTSpvtkFDa1XOiyNGHu2SPH3D47kOywZjho0PyUY5QgtZaQtiJKdKMlShBrd2ClaJOne1Pcj44mPCHDw55uDumyEvyaYFUkrgVsjvMuNWPSQLJrW5EaezK3ir9m374kMN3HzJ8uMtkP8WUBhVJkrWE7tGINSUxRWbLVpeUwnZctkdTHjtdH+2OmU5c2loKkk5UB/e1VkAo7R42UWZ1J+3w3YdMD4YMHx0wHUwBiDohUTciH01YUxLZ30Cs3MFoTaGpHf43Phzw3tMRWVrUNo5bIbtrGRvdiEhJ1pOQzXZQp+z00QHF9kNbrg93GO+M0a5HFfdjundXWSk1ydYYubYF/TuNoGrL9g8fHPLB7ph0lKFLjRBW191hRqENK0mIFIKNlkKYEjMeUOw+Inu67Wy8U7fZsBMStgK6gxFCSWRvFbV21851SxhOSx4Np3zzyRG7w4yHTq4xhjAOeNRPWO9GKBfMN9sBoGzZDg9rGx+9/4TpIK1tHCQB3bsDyrygvTVGrmwgV+5YX+E6pY8qG2+PGA2nlC4TEcYBW2stskKzEgeESnCrc7M74fwAbjHEVa59+OR3f7f5j3/6p2/wcTwej+fl5l+89RZf/d3fvbaQ9Kkf/LT5tV//zev6OgDurna/aox581q/dAm5Urej2+3ycz/3c3bpfzZGvvdlsj/9ff7srX/D499/wv7TMZPS0FKCdhyw8pE+t37gDnc+8720P/1jsPE6B/03+OOdCV/69j6/8Y2nHBxM2H8yZHw4pMgmBFGLzlqPdjfm/usr/JXvu8WPvL7K92y0SKaHiHe/zNHbX2L77T/hyR88Id2fsn+QkmlDN5Cs3G6z9X1bdO+ucvsz30fy6b9Cfv/P896R5uFgym/+2R7/9ps7bH9wyNHemDKbUBYZYdKlt9Fl87U+n/3EJh/bbPMjr6/xibWI6MP/j/R3f53xB4958tU/5enXd9j/cMjQ9ehWQsnKvR5b37/J7R/6JP0f/hHUnY+Q3vlz/OnelN/64IBf/8Y22/sTdh8fcbQ7JE+HqCBCRS1WNjts3uvzY5/c4i99ZI0/f7vDiioIP/gak6/+vzz+8h+z/bWHTPZTDp+MGRaalhKsrLdorSVsfd8Wtz/zvXTf/Mvoj7/Jg6LN01HOv/3A2vj99w8Y7IzJ00lt4/ZKl427PdZWW3zue7b4S2+s8b2bLfoHf0bxh/+G9OH7PPnyH7P9B084+PYhg6ykNNBSgrWtNhuf3ODOm99F994mrU9/jvKjn+ZbY8WXHx7ym9/c4ds7Y548OGSwNyEfDwBQUYvuWof1Oz3e/MQGP/bvbPLGSsInNxJaT77O9Pd/g6df/hrbX/uA8c6Yg/cHHOYaJWCtG9FaS9j45Aa33/wEq5/5i8hPfpan0S320pLf+uCAf/0nT3nn2/vsPRmRpznZ+BAVRCT9PqtbbbqrLX7kE5v86MfW+Qu3O2xl2+g/+R2KR99i+8t/xJOvPWDvnX32hxmlgUgK1lYT1j+xxq1PvU733ha9N/8SfNcP81is8tVHR/zmu1bX9z845ODpiHRwANhVpu2VFdbvdvmej67xVz65xb1ewg/cbrN69AHlH/8W+19+m8dvv8Nkd8z+e4ccuNFfP1LE/ZjN71nn9g9+jPU3f5DwB/4yZf8OH2Qxb394xL/+5lP+4L09dh8dkU9LpsNDAOLuCv3NNp1+zKc+vs7nvmuTH3qtz/1gjDraJv+Df8PO73yV7d/7Frvf3Odwb0Lm5pZWWyFrH1/h1qfu0bmzwdpn3kR934+y17nHH22Pef9wwm98c4c/fm+P/ScjxoeHlG5la9JfZXWrw8c+ssqPf88t7vVjfvBOjzvmAP70tzn63d/h6P0nbH/tAbvf3OPwKKvbbKsTsvbxVW5/6j63/uIPkHz6xyl7t3kabvK1JyN+491dfuebO4wGKQdPx6SDAWWREbVX6K23afciPvmxdf7d797iM/dX+K6VADV4hP7Gb3Pw9ld48vY3efr1p0wHWV22K6Fk9Y0+t75/i869TTZ/+FPEn/oco63v4Ru7E94/TPnNd3b4vXf32Ht8xHB/RJnZlZ9hu8/KZpu791f49//cbT621uaz93qsJzbN/y/eeuvZvPIcAr/IZFGeaVxtpin50ZDhkxE722M+TAuGhaYbSLq5Rn44pHd3SDYY08pSuzwZu3x8dzhlMEgZD6YMd/cY735IkU0Iky7wOrowHKwmHIxzpoW2Cy+MxkxTssGY0faQ4aMRu+OcJ9OCSWlYCSU8GdPeGKEiSTYYEbt9PoWGcV6yN5wyOkwZHqSMnn5wLMAJ+VFa3Yi94ZSNbkReLaPMpujxmOxoxHhnwuDRkA8meR3gBoVCPhrS3myTHdn9RdWG0lxrBmnOztGU0WDK0e6Q0dMPjgU4FXyUdj9mbzh1ujr7ZqlNwe4eMXwy5nBvwgeTog5wr2nDWlrQue3kTid2o70xpKXmYJyz53Qd7nxIMbGBNYhaGP06SSdCSMHuMKv3IaIL9PiIfDBmvDPi6NGQB+OcPRfguoGEp2OStYR095Co3yZxNi6NnafZHkwZDafOxg+YDveQUiGDCPgoSdums6otE5Vck45cus7WqQ8mBYd5iRKCSWlYnxS0NlpMD44w6ag+eSIvDUfTgqeDqdP1CcVkRDY+REpFkd0mil9DSlHbuJofMumIcjwm3RswfDTkyWDKk6ndUhFJQaYN0aOQ/v0BUa+DSUduFaWbkxpm7B1NGR5MGO48Ynq4g9YlKogwuiRqBewNM46mJWlb1wt5zGTE9OCI0faI0ZMxD4cZT12acb3Q9LOS+MOI/hsDK9NtPykNtdzRYMpo/5BsPCA9fGrrW7qFDO4D1DYujUEYY5fpjwfWxk9GPN2xbTbTdrP0qDQEH47o3z8iSCLMZOQWmeBWxmqeDqp6/JR0/0kd4IpsQhDeY2e9xSDNWW+FaFybzVKmB0OmB0NG2yOeHk75MC0ojQ1w/bwk2h6R7g0oh0PIpvUex7rNDtK6bNP9J2hdkrcHCPkRjDE8HUw5yorZliKtZzZ+MuLwyZi9rGTH7ZEcFIrg0Yju7Q4qGWDScb3Ke1poJs7G48GUo90Bo6cfoIsMrUvi7joqeINBP+ZgnDPulvUWm9Puz7sO/CKTxXi2bQJFTplmTAdT9vOSpy7QTErNpJSsDTOmg8xOnudZfRpDrg2TrCSbFkwnOdOjfdLDp+gio5xOCFpdhFRk0y7DNCevVl0YuxcqH02YDqYcTHIntyRzWw86SrI+mDIdRJRpZjdHQ70YYZyVTNOC6dDKLLNJLTfs9MkmHY7Swm4odXKNLsnHKfkoZbI/YS8r2ctKDnPrnDNt6AeS1cHUTmIXWb3fJq82xVa6DvdID59SpENkEBG0ukyHa2STDpOsrJeZA5gip5xY+04OUvZzzV5WMiw0kRS0lESNctYOM4qR3T4gnBPMS21tNy3IxqmVORnagB61CJIu08kqcStgkhXkWruFFwWmyClcuY5ccHs6tQF7Utrguj6Yko/dqjRd2pWFRpIWmuG0YDopyMZD0sFTsqM9hFSoqGX13egyTIt6H5U2dq9hmU7JBhPSwymDzNanagRXzd+v76eUk2o/oK4PdZ5kJdNK7tE+2eiQfHyIcIE1HW0RxgFHadEI5hqT5+SjlPRwymQ/ZcfJrUZwkRRsHU7JR1PysbMxdmuDtbGd/0pHGdPDHdLDp3bpeRAhg4hsbZ3htCAtylmHyTn9fGTnhEaDaV2nqqpeGtgaZBQjuw+xuUcrLazcqu1Mh3tMXYADSPqbpKG0nyl1Y/9d4TqlY5v5yHXddir73j6aMh1MiVftHtNKbu72v42mBek4J2u0WSEVQiqmq7fJJjlHqa1PxjR1nTAdpLWN97JZm8204db+lGzofIU7ks4YG2yO0oJsYv81bazzjLDTRwUbjCc5k6ysV25aX5FTOBvvZWXto6o9bmuTnLXBlGSUotPZqS62/RhrvzSvfYUu7LMZXRL31simHeef9NxB5teM8ItMFuXZRnC6xJQanWsX2AypK2kl3LFUWYnOqzuvXCPVhqzQ6MJQ5HoWZKqRwNT+rkt94noMtMZoTZnZ47aqgJobUKUm1YYiLSgzjSnt6Q3gNnrXcrWTN3Gyc2QQUU5TSiczc/teapmlRucFOreymrpG7ucyK9FZAaVtBJXcrND1Yhad505uitElwvV+y3K28Xyma4l2uk6qE1icjqWpOhKCMi/tUu/GaSP1OYmls2+e1XIBymxCkWuKTDtdG+WaZ+g8R+fliXIFXLk6m2QFpnQb3qsN/IWmLLXVa2pl1se05RmlW81ZbagFMKXVocxKirSo9U21QWLt3Q0MZV5S5kXtBKFxPmKpKbLCys1mcotsQlmU9rn03LmcdR22HZpK3+pot0lZ1SdXtpVct3+vLtsio3R1ymhtHXCRUeQzudWG+qrt6Lxo1GNdy7Vpflufynxm32qflnb7Dau2U9m4aju2TttFEOXc+nVT2npY5qWTaeuztbF9hjKzbaxZn0rTbDv6WNupbZylMxs3221pbVxmJWV2XG5UGiJprJ3y8libna/HRXbcxjKIZjLLysYzuUaXlFlBmZek+njZNnUtM11vRzHOT1T6FllZ+6NKrs4zCtd+qnK9aXx8W4wrBzhvaI/Hc92UzyFIvKy4oz09C/BsG72lskt7Q0lLCVouz1H9HEmBihQydKd2u/y0dBucZSAIQomKWm5+xi5EULH9XSrpLlFsICVCSlQkXZrOpeq0oaUkiRQESYCKJELZKzHAHusj3V4sGUikm/+q9JBBhIoTlJMZBXK2CVlKu3E9DAiSgFBYHTM90zepdI0CUDZlU8mNArsdQCqJDENU1KrTWCpuEUQtlLKntxy7q0wqpNO1smdLSUozS1G2lESFCqnksZPupbC6KuXsG870VVHLpgtDSRBJp2ujXMMIGYbIUNVlmciZrrZcrX1lFCCUu/RG2A3OKpAoJa1ecQuVTeoUpQwjVGBlNotWKKuDihRBEtDKSiLp7OrsHQqBChUqdOeMitlpN4HTNYgCgqiFzrM6fRZELVSg7HNJ4a5jmtnY1mFFqEStb5WibKmqPrmybciVglnZBlE9p9pMUQbhTK49C3PWdmQYNOqxpKV0vZCnajsqtPa17ad5nJys246KW6jIHsxsf24hA1fv5qaEhLL1MEiCWf1xKUr7uytbebw+KVG1WfuvqkNVejKI7DPUNm4WrrI2VpFCReqY3ErXOLJbi5ptFuyxWKou2+M2Vk2ZqrJx4yxXqezBBKEikcfLtqmriiQqnLlC6epIFEiCSNX+COyoUIYRgWs/VbneNH4ObjEWD3BC2pMlkoi4H7MWKiaxmS0yCSTdbkTcj+qTGaobdEN38kEUB8StkLi3ZlNmbpFJ3F0j6bSI4oBuEhLW1yBLRBARdlrE/ZjVlt0fV+Xxu4Fk3a1Ai/sxKonsuZRYxxu6DalxEhB318hXtijSIVqXTu46USuglwS0I1XLFYHVIewkxP2Y9UgxcPNgACuhYi22csNOyx5a6864C5Ww31Xp2l2nWBmRx626scbdFaJWQCuyJ21UjVQEIapl7Zv0YtZ0yjBSteNdCyX9Tki8EhF0WvbYI2GDVaiktV0cELUTkpUtiqhFng6trr014pZ9vxXZPVpSCHtyfRASuHLttEPWs9ncRDeQbDobh+2EIImsI3QHJyeBpBsHpK2AqN0l6W9ZZ+46EUl/kzgJ6SYBoZJOX+u8VRIT9VskKzH9tGArDoikXWSyFirWI0WylqBaESJOaplS2PoUxwFxKyDqrdkAE0ZIqYh6aySdkNiVbXUMHFIiwpCwk5CsxLTWEjbTao7GLjJZCxXJSkzYiQnbzsa4mwKUpJsE7MUBSSciXdkEqBeZJCubRK3AHp0VqPpQYIREREldnzqTknW3+AFgPVL0A2nbTiexp9U0gmoSWLlV22mSrGzZ8u5E9jNKNoJ5gIgTwl6buB+xtiOZxKpeZLIeKZJe7OpxUtcnXJtNAkUnDkjaIVFvrU7LVjaOWwFRK6SX2PokRFPXFnE/sTZ2KyjrRSZO16jrfIVUGGH/PlC2zUatgGgaEK9s2vSuLonaK85XhLRbIa1I2cDqruARQUjQaLNVHS6NYSVUrLbCWleZJPZcTCFc+xHWfklI3F2jWBkeW2QSd3tEVdtR8kaD3IsYwQkhPg/8I+xtZf/YGPMP5t6Pgf8F+CFgF/gpY8y33Ht/B/gb2Ive/wtjzK89x0c/xrON4OKEsNele7vD5q02am6bQO+1Lq3NLlHfHSnljq9KAslGN6bfT9DaMN1YRwbRiW0Cq72Y1XZIHEjbgxESESdE/TadW126d0cE+1PU3DaBzu0OrY0+Ub9TO4dA2iOp1rsxrV5MdzXB6NePraLsrrXo9BPWuzF954CtlQJku0282qO92aJ/t8tHHw3rRSYroaR3t0vndpuo10YkbXv1CLYH2k9CNnsxutCkG3aVaHMVZXc1od2LWe/GTldn3ygh6rVJNnr0XhvYDfRum0AkBesbdptAe6Nj5cYttBB2462yp5WsryQMVxN08dqJbQKdfsxqP2GjGzWcfoBs9wj7bdqbHXp3u9zPSvqBPLZNoHOrQ7KxYuVWowd3xuWtvt3M3V1N0Po+Uadv33fbBNr9mI1uRDtU9TmCyADR7hGvdune7SJDyevvD+i70WW1TaC9actBJJ1Zh0kJenHAVj/mcDUhm94mT3PCcb/eJtDux7Rd2cbViFVIRNJBtdsk6326dwfcHuVEc9sErK59W4+d3ErXjW7EMC0YrrYo8rsEjaxAe2WF7qrdiNyLlTsf0h0w3OrYE2FudZBKcC8t6LgOVbVNoPdal9ZGH5F07NFvrvNSyf2wH5OurRAkndrGcdfK7DRsrITAuM6LbPetjW932BpkRHPbBHqvdUg2etbGrY5z+lZmO5Rs9WN2+jHZxhYqiI5tE+iuttjsxfSTkFYo7bVTLsDFq12ywYjOrQ5bh9N6hWq1TaB7t0uy3kd1uxDFGGFt1Q7tCUOdvq1jRX6bIGrV2wRs24nY6sf0omAWCKSc2fh2h5XdMckgo9XYJtC926G92SbZ6NvbMZx/igNJy9m43Y/pbfSxK7vtwrGw3ae7mtB3bacdqtlddjfE84xvQggF/ALwE8AD4CtCiLeMMV9vfOxvAPvGmO8SQnwR+IfATwkhvg/4IvD9wGvA/yOE+G5jTMmCCCFK4L8F/o5xG7eFEL9rjPn0RX+7eICTEpl0UP11uvc2WBtMifsxRVoQJDaV13+9R/feFsHqOiLpoKWqD/rd6ER8ZLNNFEiMNsStkLJcsSeZ9GOiVshHNtusxAGtYJaikUmHeHOdzr1N1kZ2pV+yM0bn2p1O0KV3b53WrTVkf6M+NTyQgm4UcKsX8/pGmyIrCSKFLjqU7iSTTj/m3kabu6sJK3E4u4dORaiVDZJpSu+NLcqsJOyE9BsnmfTudu2esFtryE6/vlYlDgQrScD99RZKilpWPu3atEuo6K4mvLHRZqvbaKRCQhSjVjbsnqSjMcn+hKgTsZYWyNCeZBL3Yzr3NonWV+srbEJpR43rrZD76y0GkxylJEXeoshLglDR6kbc3uxwqx+z0XaNtHK+nT7R2hrde5vkoykqkvSGOUYbgiSgc7vNykc2aN1aRa2uI1sdCnflTj8OuLNqnfx0WhCEimka255xIK2NtzrcXW3RjVQdWI2KUL01Onc2WDkY2p51K6A/yOypIGu2J95/Y5P21qq1sQzqg35X4oD76232htYJZdOSfGpTv0knZH2jw1Y/5lY/pheperQqO31Uf5323Q1WD4YAtPdTjDY29b7WYvXj63TubJBsWV21DFACupHirtP1cJIjJKTdCGMMSkla3Yhbmx3ur7fpx4E7rgtQthPRurXKykc2iLohQko61UkmXTuy6N1bs/Wpt2rv+hOSQNnzPO+vt9keTNHa2H1wExsA4lZIb73FWieyNo4Dm6Z0F3zKTo/O3Q36uwN0aYj7EaZxksnqR1fo3Nmgc3fdHhGmgjqoriT2yLOnm/bM06QTUtYnmYT019vcX2+xkgT0IidXS2SnT7Kxgs4KVmobT9C5Pckk6kT07/dp392wR6Epe+9e1Ym4u5pwb6PNU3e+5rgd1ieZdFcSuk7XfhLUbbYq2/bWKr1765RZyXQwrcs27kd07/bo3Nukc2cD1VvDqAgpoBUoVmJ79NijjXZ90lF1kkkUB/TWW3xks816K6Tr6tPNIW74+0/wGeAdY8x7AEKIXwa+ADQD3BeAv+d+/hXgvxc2j/oF4JeNMVPgz4QQ77jv+61neJ4/wq6D+r+FED9ljNnjkjH/2c6XiWJE3HI90S5CiWNHdcWrXaJeG5m0baBxaYvQnca90o4otWE4LRBSnDiqa8U5XiWpnb6IEmTSIep3aG10UZFEKlEf1dXe7BCv9Yj6bXvqveuVSVxqJ7I9s72ePZpq/qiuqkfWCmdpByMDK7fTt2cTbnbRpTl2VFdrszsbWQRhLXfW048pteFwbANFEav6qK7VnpXbi63Dr0tOBrYXutYlXu3WKxHnj+pKVnv2FPi4RYmd52n2fqvR4/xRXRvdiPVuTDtU9f1htY3bthefbPRcQD9+VFe82iXsdZGulw+z9NlqOyQrNKu9GK0NoTvnsjqqa6Mb0U8CkmPznNbGUb9d62rl2o5CdVRXvNpDdXv2BgeXxqoc8Go7ZL0bMZkWhHFJPlX1UV3r3cjaOApmo2Snq2h13Ei5TzfNUZE6dlRXvOqyEMlstCpcfeonAeN2xFYvpshKlBv1V0d1rXcjVtuhvbfMpYFtJqJF2LPlCrj6a+tMdVRXstGfjZCrFCV2dL7SCtnoRhxOYrKwIHAHksctmy3Y6EastOzZm8rJNO6OxqjXJl7r0h2nSCWOHdVV6arabUSc1PPm1Q0fts5EjBpttjqqa6sX1222rsduKkO0e4T9CclGnyLNUJE8dlRXrWtrNloVAhLnKzbcsXlZXtY3d4dxwEovrsvW6uraTl2P28RrPdpH43oOEOxRXa0N22arrEs9pyursnU2Hse2jBpHda03dE0aWZcboTF3+5y4B3zQ+P0B8NmzPmOMKYQQh8CGe/235/723jM+T2GM+S+FED8F/KYQ4j9ldv/sufijujwej+caue6jun7o0582X/rSl67r6wBotdvfBnYaL/2iMeYXAYQQPwl83hjzN93v/wnwWWPMz1YfFkL8ofvMA/f7u9gg+PeA3zbG/G/u9X8C/F/GmF9Z9FmFEL9njMynEuIAACAASURBVPlB9/OfA/534A1jzOpFf3ulEVyn2+Vv/9zP2V90gTp4iJockH3rj8+4LmcVtXaLYOseeuUOJu6SJmvsTgqeju1hpuddl7PZjrjXT9hoBWy0FEGZog4fw+FjyqcP7bUqWXpMrup2UWu3kEmH4M4blL1blP277KX2JIadccGj4fTYdTmlNrXcjXbE3W5ML7Yy1xOFOnqCOnpSH97avC5HKLsCS3V7qLVbqK17sHoHHXcoe7fZmZTspyUPBum51+VsOrlbnZDNliKWBjV4hDx8Qvn04bnX5ai1LdTWPczqXXT/NgMTMco1+5OSh4P03OtyupHibi9mox2wngQk6T5q8AjSEcWT98++LmdlA7VxF9ldRa/cpuzd5iCX7KUFj4+yc6/LWUtC7vZibnUiWoFko6UIR09RR08otx9Q7j4++7ocV59Yu0vZv8M46DDONXuT0h7yfMZ1OStJYFPUnYjbnZD1VkCnHKEGj5HZiMLJPe26HNnfQG3cQXb6sPYaZe8WRyJhd2JvTjjvupyVJOBuN+ZOz+q6niha2QB59ASz84By99G51+WozbvIzXuU/bvouMNhGbA3KXkyynh8ND39uhx3hcytbsztTsR6S7GiCmR6hDp6gt55SLnz6IzrclYJNu4iOj3E+mvo3m0mUZ+dScmk0Dw6mp59XU4rZKMdcb+fkASC9UTRMylq8BgOHtvrcnYfn3pdjuyuojbuoG6/Ttm7jYk6jFSHvUnBk5E91b95XQ7Mrrhqh4pbnYg73bhusyIdoI62Mbu2XE+/LmcVtbKBaPdQt+5T9m6Td7bYdbpujzIeHU3PvC5nox1yv5+4A+GDehX5dR/VBZy44Pca2DnnLMqHwOuN3++71077zAMhRACsYBebXOZvr8rfrH4wxvyhEOLHsKnQC3m2jd5hgi5aqLVbdqXcZEQ4f3libxVaPUwQ2ZuvRTUfpthsR4zzkjiQjPPgxIWn64lN7VTrEIyQmDBGtrrIFXtHl8lSZDclbFx4KntriChBh636XjYl7dLqbqzY1PaSwl4UHLvwtBXaK1z6cUCkxLH0mQnbiC71fWSy3SNoXHgqWh3kygai00eHsZ1HwK4+i5VgPQlpu/vIVuLwxIWnK4mdo7GXRQrA2Hm8uI1c2bB3dGUpqjs5fglolNirPtp9a2MhEaZa9SZZa9lrU9qhOuXCU1WnY+vVjFJZXcGe2C8VstNDzV14KnvWOZggni2AkHZhxkpiT1iv0rNNXeNAshIHrCQhkRJ12VY2VisbtqyzFNkeWBs3LjyVvTVEdwUdxCAkEluf4kCw5lbVtkI1d+GpdBeeKle2s0UmJojRRtv6BLYc+6dceLqyAXEL/f+3d34h81zlHf8+Z2Z2931///wlsTEoLUID3lVosFeFUlP1oqAXKoLI78JQBHtV2qoVKtgKVqHep1WbFoRqe2EQMSQ/K3jVmtZQ7EWJtIQaopJEWyvNn9/O04vzZ87MzszOzuzsnn3f7weW992d3Xnmeeac5885Z2byws4/wz5M9MaqCE94btq40jXHwthLGUTE2ji386so1zDnL6E8u4Ks8cBTOb8Gc/1uaOEWLblFJovMzjmWpbqHqtoHngLAMjPhgac3QjuObJyv7D7Xa8jyDKbtgac37g4PPPUrCxeZQFVwl7PxMjd4+azYeODpXasCy8x+P7MNyvbZK9dhjAn6eRvXHnh6w+qKaIgyz+y8rn3w6RpnRRbk+j67co+ZWubR5Qkmtza+dtM+SSEvbP9pPPDUXHsNzOqKlZv5+XqEc/vqWu2jnZbFxgNPb57ZPr3IZPZ7Rcr+A1wf3wFwv4i8ETY4vQ9Ac+juUQC3YOfW3g3gm6qqIvIogC+JyJ/DLjK5H8A/jTkIEflDVf2Mqv6ziLxHVb8CAKr63yLypiH7mHaZgGtE5upr7Pza+bWqs2R2vF/Or0LzVZg4BqzzXbgFCavc2OXW/qnWImG+4uoyqzqKk6neEV67aZcT+8fQ1xzSdaBYWCfoGm0m9YZrRMJTxAEEuVcXtqMURqpriJyeAGzGlxcoVz+v7vTgnf6V6yjzFTRbhHmpzMm8usywXFsH9VJRPQnZO/1ryxyr2OE7uWVxZp0D7EMoa0+bLha2o169AXVybaAR5GrlXl/ayfpVZkIH9brazmmwzKK5MDHQvABwDrl6B5nJ7BO93a3HYhuXxZm1b7CxTSLOC7il01WA27Bx5q5diwNrVkBWV2DKtX0Q5dmVIFeWZzawnF+3Dj/3c7qCTNTa2M2vLPMs3AbMvz8vDFZ55pygu67IGLsfLW3S5OToy/9XNfPlmZ2fWV2x5zWzSUQGa99lJjBujnGVmZquRWYDq7ezTyLsvNTCntvrd0NfecnOd3ldiwWQFzCrK8Dy3J6PcJmAhOTFJxKvumBu+5Y9r0VmvxP6j5uD02IJWdpEoiwWVZ91cmV1xSZp2QLIF9G1hjbYnBUGN9SuMI5t7OX6Plu4OV0VY22WryDnZtPGvh2vziFnV1Fmbv7a2O8WLnm5sczxamGwzDO8fGcdbOF1vbrIGsHcntuyOAv+Cd5XeF3dvLpmhX25PuD7zzKzSVORCV66U9d1mZuQCOchKb0YuDm13wXwGOxlAl9Q1X8TkU8CeFJVHwXweQB/4xaRvAgbBOG+92XYBSl3AHx4wgrK9wH4jPv/YwC+Em17B4A/2raDaYtMTG47wfIcWZZBV6+E20/BZEBhG3booMZOPFunb1CiRGZM6Cxrtdfj+OtKFj4TjNtOtoAWS5TlHchVgVFF6TqLuM6i+dI6+3xRuzjWBxtVQFxw9dceeeecGwnVm3f6ajKrJ4BSS8hVQbZYhdtxicmAxdIGt3wZMlDALkbIjWDl92kKnK/LDbmLTJBnlcMPiwLyBXS9glwBTLkONgYAszyzy7/zFbRY1hbUZFI5BxE76V+q1mzsg7j/6xcF+BWgKM8h57ZCjXUNNs6X1sYumEuURBixNn4lr+vq5XonWK8uFlAnE+UdyCuRjfMFkOf2/DsnqA2nv3AjA4WpdAW8Q3IX7zqHBFjnK64Na2mrVskX0NV5aG6SL4DFyiUQRVVdoNpfJgDcAqFXo9tMrfIs2CM3iBImZ+N8CVmWwdn6+6YGG2dFLXHxFZzfJ2CvrbP3Wa0CnK+eV7k9Nr94CHGwWZXIigX01VfqNvZ9NiusXHetYWYEBWyFiIWVE+vq5a584iLe6VeJMMS49lTZuOYrinP7UORwHZyE/qiFwaulWLlFdYH/IqtGDuJrSP25VWfjLMug63XNxr7PIlu4hMmuYPZ9Z5EJSjWdfdYnwmG19WzoHEOU/RJVvw7g643P/jj6/yUA7+n47acAfGoPhyEd/7e9b2XaEKWIbbzZwt2fzwYcH1RsJpbbqiu+wwDcXS/E9r5S7TBA4bZZZ2s7ib+Q19hRO6gIxFWPyBbQ8g7M8iwcj3p5JgtDZ0C1sjBzGb/tpyYsxfFyC2NXqtWCqnMOMCbohAVg/E2RnR3sMeXV9+ECq6irQtXe8aMhN3eO3usbyw1O2Dt0Z2PABV4vV4w9DjFhfVFmgKy0+1br5Zx8t804e5iotRgTKmVxNm7q6m3sj69a4ecTCXuvRxtINm3s5Yo7twHfTkxuz3Es1w3RabaobAyEu4PYfQsyBYoMUJUgtzrv9twa1xasMayuMBmQLexFtVl1Bw/fjtRkNbm+Tfk2XASbx3eTqexrRIKNrfM2QJZD17n9PC+DXBWpKh9/55To7h7expk6FSBBriA6py4wVT/059bqCnXBNdK12Y41snPm5ZaCItvUtXAVedZc9ed11dKe22JR0xVNXSOZRnx/1HC1c03X0IYl+JUNuescyEo70hLLdaNKQe/IviIa9O3us5Wd549xl/JeZtrxf9v7VqZVcH6YMqueuK3NbX7J/EbGbc9ZBusY43vR2Q4iUfkf7dM7ucI9mkJLqHvMRbW9sB00CnC2qvHDWQYiirWpbrrr/BTykG3Xhztg8rpu5R3bYRu6BmdoqsDqhzy8HqVBTW4mVSaYRY7QdtAFtKjuQl+zsQ8yeVFl+bC2s8ENcKN2yEy9j/ghJwMJDgKIqhrADVVaORpnkO6c+6w33DlCgMwoCrX7FZdIdNo4doTexsXK6ap1udG5rTlDf26djb1uJSq5mfjz6eXGQ97u2rTcPlkcYjZ1dRV8OLdu7s9XU6IS9Gva2N9BZyPIeaeaL4Ey79bV2Vhr59ba2Ove1p58ZW5cXwI0qmp8QmqqvuN1z4q6rqiCJpxczaxucZ8VF+DyKMiFpNTpKrnzSiav91kgyPXDk/UREKCETV5iG3td/bx+ZqK5sDiQ55Gv6OizQW7Yr6uUYfe5jgJY3Gd94jQrevgKLhF+RUT+B7YZnLn/4d6vhuxgpwBXO49xVubmQ7DONrfHjScaz88MkCugsM4wvou+z+59JzFNub5qBOydx30m6TJdL7MeaAQCdbffsjlZhqZDcoEhZKGV04dzhFZfAGVWNTov189HxZkvqiEPozYbLdEMrNUQVi0DNaZyhAWsvKaNvVxTr2pU7f5KNxe31nZdfRZs3O9ipx+C5jqrdzB/Xv32qLrwHV7MdhsbqRySr07DkKGW7TYO7ck2XXHnNhMNQ0Xehm2OMHMJU1WxusTFVW/qHGHNxlm+YWNrZxe4Ijna4vT9cHczUUO2qJ4jtu7WtRbMnS7hBjAt7anm9BujAmoyG+SKlX3+WV+fjaopEUC0ki0CZA1dK6cfVY6xrmKnKHR9p95ngWBjPzwJ8dMZriqHbcvNAGcTKQS5tbkwkwNmHeZrUZbtfdaNfsR9Vn2S4Gwc+6dmn40Tl7k48CKTJFDVbPu3+pk4ROmHWhbV/9G2EIx8owVCRplJlHkaDVmwiZyPn0uoN1oD1QyizjGZMpz8WGbIghtO33caO1wo4VEXXm6VfcYpW1XBCRYu4LXLtUNssWPwjtXqKD1yBfVgXrNv0D9vlxsFcy/X2xilBn1KrWSakEg0bCwGalDZWMymTJ8hN2wsLnmx57NdVyN+SDZySP68qVZOv8vGXtfY6RupBoCdvFhuLWmqnVdjx2/Npq5BbhiBkCphgqsUBbBN18qNH50SD1GKYCNRU+Pasaum2mxcC6piAK0qDBWBmH4bbwRVv0+4NK/ZZ2Nd46Fgt/PYxl26+sSlSiIMtLTJS2t7AqqkRWpWCgkijJ07j+XWdA1JU/3c2mStv8827eyT4djGWaevOEAFB1zWCm4y04YogVDFhew3xneUuIM6qk5qnWEz8/WNqDa+7ZypxENLQH1IyTted2x1mW5lpPPyqqg1zqpqjIZYgi4GgA9ygH+a9IbcaDgJ8B3fiow7apvczK3GqnVSMVApIX6eM5brh3fiQCP2MoEQbAT2Ti/ORF5u08axHjbbRWXjFl03Xl5fuGCCbhv7edWars7GqpXMTl0bNo7l+mQplhvr6m3sUTEQ34b9kG+zgmsJ5na/9vhL1RDU29px5h1+0LceWAHUh56bchtDdlWFsakrgLCYZSPQeKfvbdylq+uzzaTUuAVKsY2buvrAv2FjU/mJXhv39Nk+X+FtHFP5CvT7isYohE+G+2zs+5Y0dJ2HSztEOZnpc3Bxyd/1HdNwglEDssvW7ZxJvN1nZht7jbNuaZlnbAmmTYzYTLTZZMKoSl97lR65W2S26erldt5rzumqKLFhjVhXsynfO0LvgJsZd2zj2jCWGLurHhuHxCWuGv1fsUM4XTbeCKoNXVR6dN2CyyE621P3D7ttrM3fuUoKqHTta8feCW7Y2Olq99MhtxFU7XdtAuETmTbXF8uNieV26trTlv1Q5WAbN/ajIsNs3KFTl9zuH5kw39i6Lda1EVj7bOz77LBWOREFA9xIJldwcXbmxmoqQka42QyqYQC3bD/KX8PjYvqyI5eZuYPYkOmPrZb5usBaCpBBajK9XD+UJFJf7BGcfmn9fpeuIdBEcvt0bcptrsaKh37tqF+H3DYTObneAccZqJfpv9ckzn77dG1zvkZ1kI1bdXWBtVfXyMa+IBcRuKfSbejaJjdUUmJqQ3WKEpDG0L8YW1n4/8M+q/bUZWMnZZOQHOS2Um1+Kda1BeMOYFt7av+xldvXZ9uS0hLVKuk+G29U5kEPA6zvdNs4GoXw51bhpzGsrl029lWjNPdZ07Uht6cd+3O7zcZe11rygn2jtaeck+FMr+CAMKTVuS1+jyrj8p94BxG2S317s9HGwWaozGrflQM2jdwvDm4by369I/Ry26rXlo4CRJ3F6eoDXVPX1g6qZbt92+RGQ0phvqBFZqxreN8nN9bVb4917XCEfXJ9UN04twAky22H3tHGJWyQC4uDGqs3Ecmt4R3hEF0jeWttBPQeXWuVhpNRa8d9urbI7rJxsz1tbcdNO8Ryo8+q6qmycVNXL9fLbFasQHRu20Z+WkYhgMrGcAFlkK4TbCzw85PtvqLpnw7BZVxksg+mz8Ehyrxh3KpGU9vWNgTgOwuA0HhLNJxuWwd1+wmdFECzqbU5QKDeSb1MoC63U+ZQuR0VlXesJvTOze2Dde2zcURN1w6Z/nuxvLCEvEdul+NdB+dTyS2bwbTp8BtU7cnJbWzbOhwcBG3quyG36Qi7dPXfbTjCoTbeIFSMA9txS/IS27hP19qcbrTIIvTZBttsHAectm0bh+P3U7NzRztu7g+VjVXVXbTeLnPzx+Ns3JTbp+vWIdJ9wQA3iukBLmpE/m1Xo2niO2vs/OPG2tt4xHQPBbQ4pFpVEwVX32manaSrgw6SG/3v5XonE8tt6tq3T218NMTGcSftotXGjSAHuAyyrXpqvI8dai2BQY8jau7HOcKt34t08DaOdT2EjYP8XW0c0amrl9ujc1WtDrCx3x+ws42BqD015O7SZ2O5rTbekqx5uuT2JYnhbax7S/K98fMOGx8uuClq5TIZzF4quJjWTtORmbkh7lpG6okbatswy9iMxssN7xsOYbLcPmfVI9cf2xB5Q2y8EdBb+kenw2+R3RvYtui8k41bZLdu72Bbe+p1SI0KZ0PeQLnNz3uJk6eubQNk9tm4U+4ImX1yB7WnbTZu+ymG+4oh9CYU6O47R73fJCu4Ucy/ijL+rqPpDLoaaOvHQ2VG29pyn75OsbUZ76BrvE/fSfvodPhDGCHXf7ytutll2xBdt2X5g3XeUS66tm+T2xHk/K5G2XhIWz6Gjfsqxi1ye/vsBLnx+8Fym7L7vtOxzyHt6RBwDm4c+6ngRmaDsXPo2tYrE9hJbrOT9h1T6z7jjtKUO8Ahd8kd1H/6qpoBdt5J7jZde2TuYuNOhz9E7han1CWzdX9xkBmh6xC5+9R1m42l8d1Oudt03afcpvw+uT37bsodVJEfS9e9wevgxrKfObgRGVLTOXSxdQhrRGa2TeYkucfStWN/Kdt4KzvoO1TXTjkjbVzbzehfjpd7DJnb5A4elt112xS5x9KVHJVRAc5ebZQ4A7LuLuZotEPkTmKCvseQuXOGP7fcbfTIBPrlJqfrnuY4d5YJjJKbcns6GKzgRrH3RSYbzNFZ/H5HDtlNYoLcPn1H6zpA7rZjGi13y35ncQ4J6krmZ1J7mslXHKx6u7xPE5jMUefggPYMeHDDacsKR8rcSXabvgM7SltHHS1zoNzJNj6GrsDO5zXe/950HSh7tNwJc7pd7Zg2Hih3bhvvCQG4yGQk+6vgJlZNkxrMSNnHkDlJ7rFsfKwslzaeV+ax5F5QG886mslbdY1i/iFKQgghE+CF3mM5eIDrWz5OTheeU0JmQsE5uJHsHODoyAgh5LBwDm4cHKJMkOQvwSCEHBCuohwLAxwhl4CTuHaVdMMAN4qTCXDsnMPhPOfFhOf0kqIKlOtjH8VJcjIBjpB9MjVhYhJBDonyMoFRMMARQkjSsIIby1EC3Kllv6c4fzHGxsfU8RRtTObl1PzEbCgY4EYy400bCR02SYlDBwu2f3JsjjZEyeysnVN0CjyP80Mbz0+qNlYodM0KbgxHnYMbGuSmOP1UGy3Z5NSGKZmkkYOg4L0oR3L0RSZ0EvNB27ZzzCB6SgGcpAIXmYzl6AEO6HfEqTiEQ1QXqei6C8cMoAzeu3OoKnlfMpikAVCFMsCNIokAB5ymc98nc+k/p4OYY79zO+B97/sQDvjSO/gdGdOGkrcxhyhHkUyAOwXGdJwhDjCV4L6LfnM6BM7FzcupJRGEFdxYGOCOzKGGi/bpgA/hzOdwwnM79bkq2lNiDhvPnUAkb2NeBzcaBrgD0NVBD5npnlqVAew3yJ1iVcEqeX5Oo08ohyhHwgu9d2Rsh5CW16EZIrNLP41eh2Qf8g5l633KOQ3HW+fUhj1PxsYK6Hq919cUROQuEXlcRJ52f292fO+W+87TInIr+vxbIvLvIvKUe/3CpAPqgQHukjE0yDVfx2TsMRwjkdiHzEPZe59yUluZuU23Y7fp3XCXCezzNY2PAritqvcDuO3e1xCRuwB8AsCvAXgLgE80AuH7VfXN7vXjqQfUxagAd9mHNk6rc2xyquevGejaquJjVsgxux7DMSvkKYugDm3rXeS16ZVCwrYzmlyAeyeAR9z/jwB4V8t33g7gcVV9UVV/AuBxAO+YKnhXOAdHTo5Tmj9qHmeqztUfV5ddU7P30DnlVO29KzM8LuceEXkyev+wqj488Lf3qupz7v8fAri35TuvB/Bf0fsfuM88XxSRNYC/B/CnqjrLqbrQAe6iNO45OMVFJ2R+Ti15uBxteJY7mTyvqg90bRSRJwC8rmXTx2tHpqoisutpeL+qPisi12AD3AcA/PWO+xjEhQ5wc3JKjoCQXTiltu2P80IHuiNcJqCqD3ZtE5Efich9qvqciNwHoG0O7VkAvxG9fwOAb7l9P+v+/kxEvgQ7RzdLgOMiE0IOxCk54VM6VuB0AvIF4VEAflXkLQBfbfnOYwDeJiI33eKStwF4TERyEbkHAESkAPDbAL4314GygiOEkIRR6BxzcFP4NIAvi8gHATwD4L0AICIPAPiQqj6kqi+KyJ8A+I77zSfdZ1dgA10BIAPwBIC/mOtAGeAmsG1inszDKdqbFREZTWJ3MlHVFwC8teXzJwE8FL3/AoAvNL7zcwC/OvcxehjgCCEkafi4nLHsFOBOLQs9FM1rs06Fy7MK7bjQxvNzoW3s7mRCdudCV3DHcOBzP6F833QdS8oO49j39dyVtmNL2b7A5vGlbF/goidrvBflWHYOcKk39CYp3jZp7G2nDsk2eak5k1OrolN2yKnaL0V7tdlqFvtxiHIUOwW4+MRNuVYmxYZ6DI7tSFI4D8e2gScFWwDp2ANIxyYxKdkHOJCN+ETv0ewU4P7lu999/uz8/Jm5DoYQQi4Av7TvHSZ2mcDJsNsiE9XXznUghBBCWlCFrhngxnChF5kQQsipowoGuJEwwBFCSNIkdyeTk4EBjhBCUoYV3GgY4AghJHEY4MbBpwkQQgi5kLCCI4SQhFFVlLxV1ygY4AghJHG4yGQcDHCEEJIyvA5uNAxwhBCSOAxw42CAI4SQhFHldXBjYYAjhJDEKVnBjYIBjhBCUoYXeo+GAY4QQlKGi0xGwwBHCCEJo+BlAmNhgCOEkJRhBTcaBjhCCEkcBrhxMMARQkjKKFByiHIUDHCEEJIwCg5RjoUBjhBCUkYB5c2WR8HH5RBCCLmQsIIjhJCk4a26xsIARwghKcM7mYyGAY4QQpKGi0zGwgBHCCEJo8qbLY+FAY4QQpKGc3BjYYAjhJCU4RzcaBjgCCEkZRTQtR77KE4SBjhCCEkYhXIObiS80JsQQlJGAS11r68piMhdIvK4iDzt/t7s+N43ROSnIvK1xudvFJF/FJHvi8jfishi0gH1wABHCCGJU651r6+JfBTAbVW9H8Bt976NzwL4QMvnfwbgc6r6ywB+AuCDUw+oCwY4QghJGHWLTPb5msg7ATzi/n8EwLvaj1tvA/hZ/JmICIDfBPB3236/DzgHRwghKaOa2iKTe1X1Off/DwHcu8Nv7wbwU1W9497/AMDr93lwMQxwhBCSOHsYVmxyj4g8Gb1/WFUf9m9E5AkAr2v53cfjN6qqIpJU9I1hgCOEkMvH86r6QNdGVX2wa5uI/EhE7lPV50TkPgA/3kHuCwBeIyK5q+LeAODZHX6/E5yDI4SQlElvDu5RALfc/7cAfHWwKqoK4B8AvHvM73eFAY4QQhJGAZSl7vU1kU8D+C0ReRrAg+49ROQBEflL/yUR+TaArwB4q4j8QETe7jZ9BMDvicj3YefkPj/1gLrgECUhhKRMYotMVPUFAG9t+fxJAA9F73+94/f/AeAtsx1gBAMcIYQkDu9kMg4GOEIISRjlvShHwwBHCCEpwwA3GgY4QghJGt5seSwMcIQQkjLuZstkdxjgCCEkYRSz3MnkUsAARwghKaPKJ3qPhAGOEEISh4tMxsEARwghCaPKIcqxMMARQkjiaMkhyjEwwBFCSMroXp7CfSnhzZYJIYRcSFjBEUJIyvBOJqNhgCOEkIRRgJcJjIQBjhBCUoarKEfDAEcIIUmT1vPgTgkGOEIISRhVoFQGuDEwwBFCSOKsGeBGwQBHCCEJowA4QjkOBjhCCEkcVnDjYIAjhJCEYQU3HgY4QghJGFVWcGNhgCOEkMRhBTcOBjhCCEkYhbKCGwkDHCGEJAzn4MbDpwkQQgi5kLCCI4SQxGEFNw4GOEIISRiuohwPAxwhhCQOK7hxMMARQkjC2EUmjHBjYIAjhJCE4SrK8TDAEUJI4rCCGwcDHCGEJIxdZHLsozhNGOAIISRxWMGNgwGOEEISRgGUxz6IE4UBjhBCkob3ohwLb9VFCCEJ41dR7vM1BRG5S0QeF5Gn3d+bHd/7hoj8VES+1vj8r0TkP0XkKfd687Qj6oYBjhBCEsZfB7fP10Q+CuC2qt4P4LZ738ZnAXygY9sfqOqb3eupqQfUBQMcIYSkO4SMzwAAAV9JREFUzJ6rtz2syHwngEfc/48AeFfrYaveBvCzydImwABHCCGXj3tE5Mno9Ts7/PZeVX3O/f9DAPeOkP8pEflXEfmciCxH/H4QXGRCCCEJM9Otup5X1Qe6NorIEwBe17Lp4/EbVVUR2fXgPgYbGBcAHgbwEQCf3HEfg2CAI4SQxDn0hd6q+mDXNhH5kYjcp6rPich9AH6847599feyiHwRwO9PONReOERJCCEJk+Aik0cB3HL/3wLw1V1+7IIiRERg5+++N/WAOmUpr68ghJBkEZFvALhnz7t9XlXfMfJ47gbwZQC/COAZAO9V1RdF5AEAH1LVh9z3vg3gTQCuAngBwAdV9TER+SaA1wIQAE+53/zvZI3ajpUBjhBCyEWEQ5SEEEIuJAxwhBBCLiQMcIQQQi4kDHCEEEIuJAxwhBBCLiQMcIQQQi4kDHCEEEIuJAxwhBBCLiQMcIQQQi4k/w8lIUpjUjlBxAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# make a new simulation\n",
    "simulation = Simulation(omega, eps_r, dl, [15, 15], 'Ez')\n",
    "\n",
    "# add a waveguide mode on the left side \n",
    "simulation.add_mode(np.sqrt(eps_m), 'x', [17, Ny//2], Ny-50)\n",
    "simulation.setup_modes()\n",
    "\n",
    "# solve the fields\n",
    "simulation.solve_fields()\n",
    "simulation.plt_re(outline=True)\n",
    "\n",
    "# compute the transmission\n",
    "fld0 = simulation.fields['Ez'][20, Ny//2]\n",
    "fld1 = simulation.fields['Ez'][Nx-20, Ny//2]\n",
    "T_linear = fld1/fld0\n",
    "print('linear power transmission of {}'.format(np.square(np.abs(T_linear))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "solving 0/10 with source scaling of 10.0\n",
      "solving 1/10 with source scaling of 16.68100537200059\n",
      "solving 2/10 with source scaling of 27.825594022071243\n",
      "solving 3/10 with source scaling of 46.41588833612777\n",
      "solving 4/10 with source scaling of 77.4263682681127\n",
      "solving 5/10 with source scaling of 129.1549665014884\n",
      "solving 6/10 with source scaling of 215.44346900318823\n",
      "solving 7/10 with source scaling of 359.38136638046257\n",
      "solving 8/10 with source scaling of 599.4842503189409\n",
      "solving 9/10 with source scaling of 1000.0\n"
     ]
    }
   ],
   "source": [
    "# add nonlinearity into the system\n",
    "chi3 = 2.8*1e-18\n",
    "simulation.nonlinearity = []\n",
    "simulation.add_nl(chi3, nl_region, eps_scale=False)\n",
    "\n",
    "# modal source amplitudes to scan over\n",
    "Ns = 10\n",
    "srcval_vec = np.logspace(1, 3, Ns)\n",
    "pwr_vec = np.array([])\n",
    "T_vec = np.array([])\n",
    "\n",
    "for i, srcval in enumerate(srcval_vec):\n",
    "    print(\"solving {}/{} with source scaling of {}\".format(i, Ns, srcval))\n",
    "\n",
    "    # remove the source and make a new one with amplitude of srcval\n",
    "    del simulation.modes[0]\n",
    "    simulation.add_mode(np.sqrt(eps_m), 'x', [17, Ny//2], Ny-50, scale=srcval)\n",
    "    simulation.setup_modes()\n",
    "\n",
    "    # solve the nonlinear fields\n",
    "    simulation.solve_fields_nl(solver_nl='newton')\n",
    "\n",
    "    # compute the nonlinear transmission\n",
    "    fld0 = simulation.fields_nl['Ez'][20, Ny//2]\n",
    "    fld1 = simulation.fields_nl['Ez'][Nx-20, Ny//2]\n",
    "    T_vec = np.append(T_vec, fld1/fld0)\n",
    "    \n",
    "    # compute the power transmission\n",
    "    pwr = simulation.flux_probe('x', [Nx-20, Ny//2], Ny-50, nl=True)\n",
    "    pwr_vec = np.append(pwr_vec, pwr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASIAAADkCAYAAAA8eCmrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXd4VNXWh99FKEFaIIDUGLr0DoJIUykKAgIqclVE4WK/Fq4oiMpVwfapWK6iIpZrofeiKAoqKCCShCYIKAm9JLQEUtb3xzkJkyGTzElmMslkv88zT87ZZ5991s7MrFm7/baoKgaDwRBIigXaAIPBYDCOyGAwBBzjiAwGQ8AxjshgMAQc44gMBkPAMY7IYDAEHOOIDPmKiCwTkTs8XIsUERWR4vltlyGwGEdUSBCRJ0RkmVvaTg9pt+Svdd6jqn1V9WNflysiz4jIZy7nNUVku4hMFRHx9fPsZ4SJyHQROSgip0TkDxEZ53JdReSMiJwWkTgR+T8RCbGv7RWR8yJS2a3MTfZ9kf6wuaBiHFHhYTXQ2eWDXB0oAbR2S6tv5y2yiMhlWP+Dhar6oDqctZtVROYhSnsNKAs0BioANwC73PK0VNWywNXArcAol2t7gGEuz2gOXOLE1mDBOKLCw3osx9PKPr8KWAXscEv7U1X3A4jIGyKyT0ROishGEbnKTq8hIokiUim9cBFpLSJHRaSEfT5SRLaJyAkRWWF/udPz9hKRHSKSICLviMgPInK3fc09MsnU3BKR713yhojIK/ZzdwPXu1ZYRCqIyIcicsCOKJ5Ld7qeEJF6WE7of6r6b2/KEpERIvKTiLwmIseAZ7JKy+Jx7YHPVfWEqqap6nZVnZ2VXaq6HVgDNHNJ/hS43eX8DuCT7OoXrBhHVEhQ1fPAL0BXO6kr1gf7R7c012hoPZaTqgR8DswSkVDbUa0FBrvkvRWYrarJIjIAeBK4EahiP+cLALspMRt4AgjHcoSdc1mtUUA/oDXQDhjidn0GkIIV5bUGegF3Z1NeXaz6v6eqEx2W1RHYDVwKPJ9NmivrgOdF5E4RaZCNXYhIE6wfik1u95cXkca2U7wF+Cyr+4MeVTWvQvLC+lWeZx9vBhoAfdzS7sjm/hNYTQWwvoTf2ccC7AO62ufLgLtc7isGnAUuw/oFX+tyLf3eu11s/MzleiSgQHH7/HuXvN8BY1zy9krPi/XlPweUdrk+DFiVzf/mJBAP1HO7lm1ZwAjgb7d7LkrL4pmlsRz2RiAZq1nW1+W62jadAP4EngOK2df2AtcAE4DJ9vv4jV13BSID/XnLz5cZnShcrAbus5tUVVR1p4gcAj6205rhEhGJyGPAXUANrA93eSC9c3QO8Kbdr9QQSMOKfMByOG+IyKsuzxagpl3WvvREVVURic1lfTKVBfzlcnwZVlP0gEtfczG3/O4sBA4D34lIV1VNL8+bsrIqN7tnoaqJwAvACyJSHhiHFXVGqOpxO1sbVXXvN3LlU6z3rA5FtFkGGEdUyFiL1Sk6CvgJQFVPish+O22/qu4BsPuD/o3VSbpFVdNE5ASWQ0FVT4jI18DNWJ2tX6r9U431BXxeVf/nboDdBKnlci6u58AZMne4VsumPgeA2i7nES7H+7CimMqqmpJNGZlQ1UdEpBQXnFGcl2Vl1aHtdSe3/T68gNVkrQMcz+GW9Pv+EpE9wHVYPxpFEtNHVIiwf4E3AI9wIXoBq5/oETL3D5XD6hM5AhQXkYlYEZErn2M1tYbYx+m8CzwhIk0ho6N3qH1tCdBcRAbaHdD3kdnZ/A50FZEIEamA9cX0xEzgQRGpJSIVsSKK9LoeAL4GXhWR8iJSTETqiUi3bMpL536sjvxvReTSPJblERF5SkTai0hJEQkFHsJqGu5wWNRdQE9VPZMXewozxhEVPn4AqmI5n3TW2GmujmgFsBz4A6vJk8TFTY2FWP1MB1V1c3qiqs4DXgS+FJGTQAzQ1752FBgKvAQcA5pgOcdz9vVvgK+AKKy+k8XZ1OV9287NwG/AXLfrtwMlga1Y/SyzgerZlJduvwKjgV+BlXYHe67KyulRwEfAUWA/cC1wvaqedlSI6p+quiGPthRq5EI0bjA4R0SKAbHAcFVdFWh7DIUTExEZHCMivcWaVVwKa9RIsIaiDYZcYRyRITd0whqOPgr0Bwba/VcGQ64wTTODwRBwTERkMBgCjnFEBoMh4AT1hEYR6Q/0L1eu3KiGDRsG2hyDocixcePGo6paJad8RaKPqF27drphQ5GepmEwBAQR2aiq7XLKZ5pmBoMh4AS1IxKR/iIyLSEhIdCmGAyGbAhqR6Sqi1R1dIUKFQJtisEQHETNhNeawTNh1t+omT4pNqg7q7MjOTmZ2NhYkpKSAm2KIUgIDQ2lVq1alChRItCm+IeombDoQUi2564m7LPOAVrclKeii6wjio2NpVy5ckRGRiL+0VY3FCFUlWPHjhEbG0udOnUCbY5/+HbSBSeUTnKilZ5HRxTUTbPs+oiSkpIIDw83TsjgE0SE8PDw4I6wEzzo33lKd0BQO6Kc+oiMEzL4kqD+PO3+AcSDu6hQK+t0BwS1IyroTJ06lcaNGzN8+PBM6ZGRkRw9ejTbPO75DAa/cO40LHkUPrkBLgmHkFKZr5coDVe771PgnCLbR+SU+ZvieHnFDvbHJ1IjrDRjezdiYOuaeSrznXfeYeXKldSq5fkXxZs8BoNf2PsjzL8X4v+GK+6DnhNg+2KrTygh1oqErp6Y5/4hMI7IK+ZviuOJudEkJqcCEBefyBNzowFy7YzGjBnD7t276du3LyNGjGDFihXExcXRqVOn9B0gMuUZOXIkt99+O8OGDbson8HgU86fgZXPwK/ToFJduHMZXNbJutbiJp84HneK7BKPbdu20bhxYwCeXbSFrftPerx/09/xnE9Nuyi9ZEgxWkeEZXlPkxrlebp/02ztioyMZMOGDUyaNInKlSszceJElixZQr9+/Thy5AiVK1fOyFO5cmUefPBBj/kMBQPXz1WhZO9PsOBeOLEXOt5jRTwlc7/5rLdLPII6Ikpf9Fq/fv08lZOVE8ou3SmrV69m7lxLrvn666+nYsWKecpnMDjm/FmryfXLu1DxMhixFCKvzLfHB7UjUtVFwKJ27dqNyi5fTpHLlVO+Iy7+YgHCmmGl+eqfnfJko8EQcP5aa0VBx3dDh3/CNU9DyTL5aoIZNfOCsb0bUbpE5i3XS5cIYWzvRj4pv2vXrnz+ubWbz7Jlyzhx4kSe8hkMXpGcCCvGw0d9IS0F7lgM172U704Igjwi8hXpHdK+HjVL5+mnn2bYsGE0bdqUzp07ExERkad8BkOO7PsV5t8Dx3ZB+7vhmmehVNmAmWM6qw0GH1LgP1fJibDqeVj7NpSvBQPegrre7zPpdBqL6aw2GAyZ2bfejoJ2Qts7odd/oFQ5r2/3xzSWdIwjMhiCjaiZmScddn8Cju6An9+EcjXgtnlQr6fjYl9esSPDCaWTmJzKyyt2FD1HJCJ1gfFABVUdEmh7DIYCRVZSHQvuAxTa3AG9noPQ8o6LVdUsR44B9ntId0K+jpqJyHQROSwiMW7pfURkh4jsEpFx2ZWhqrtV9S7/WmowFFKykupAoUxluGFqrpzQkVPnuOez3zxerxFW2nGZ7jh2RCJSRkRCcs6ZJTOAPm7lhQBvA32BJsAwEWkiIs1FZLHbq2oun2swFA08SXKcOea4KFVlwe9xXPvaD3y34zD9W1QntERml+GraSw5Ns1EpBhwCzAcaA+cA0qJyFFgCfCequ7y5mGqulpEIt2SOwC7VHW3/bwvgQGqOhno52U9srJ7NDAaMMPchqJByjmr8/lcFsuVHEp1HD6ZxJPzYli57RCtI8J4eUhL6lct65fF3+BdH9EqYCXwBBCjqmkAIlIJ6AG8KCLzVPWzXNpQE9jnch4LdPSUWUTCgeeB1iLyhO2wLkJVpwHTwBq+z6VtBkPhYP/v1ojYuZMgIaAuncoOpDpUlbm/xfHsoi2cS0lj/HWNGdmlDiHFLK2lga1r+mz+nCveOKJrVDVZRMqmOyHb4OPAHGCOiOSbSK+qHgPGeJPXV2vNDIYCS8p5WPMKrH4FylSBW2dCUkKupDoOJiTx5Lxovtt+mHaXVeSlIS2oWyV/Jjnm2Eekqsn24Sa77ybDeYlIQ7c8uSEOqO1yXstOyzNFaRePsmWz/8DEx8fzzjvvZErr3LmzP03K4Ntvv+W2227zKu+YMWP46aefMqU98MADXHbZZT616b333qNatWq0atWKunXrMmPGDK/uue+++3xqR544EAXv94QfXrQczX3roGFv6/jhGHgm3vqbgxNSVWZu2Me1r/3Az38e5al+Tfjqn53yzQmBs87qSlidyrEiEi0iM4GlPrBhPdBAROqISEms/qiFPijXt/ua+WkblfwiK0f0888/58uzN2/eTOvWrb3Ku27dOq644oqM871797Jq1SrOnz/PqVOnfGZTdHQ0zzzzDL///juzZ8/m0Ucf9eqe5s2b+8yGXJOaDN9Pgfd7wJnDMOxLGPQulHauxrA/PpE7PlrPv2dH0bh6eZY/1JW7XJpi+YUTR/S3qvZQ1WpYI1zvAyOcPExEvgDWAo1EJFZE7lLVFOB+YAWwDZipqluclOsJn0VE6XMzEvYBemEbFR84o4EDB9K2bVuaNm3KtGnTAOvL17hxY0aNGkXTpk3p1asXiYmJHvO7MnHiRF5//fWM8/Hjx/PGG28wbtw4/vzzT1q1asXYsWOBC1HUJ598QosWLWjZsqXXkcuNN97IhAkT6Nq1KxEREaxcudJj3nRHdO7cOUaMGMGTTz6Zpajbtm3baNiwISEhFwZln376aSZMmECTJk3YssUnHwsAoqKiuPzyywGoVasWqampOdxh3RNwR3QwxoqCvp8MTW+Ee9dBo76Oi1FVvvj1b3q9tpr1e47z7A1N+XLUFURWzv8FrxkGefMCdgOdsSYSen1fIF9Af2Ba/fr11Z2tW7deOFn6uOr06zy/JlVRfbr8xa9JVTzfs/Txi56ZFceOHVNV1bNnz2rTpk316NGjumfPHg0JCdFNmzapqurQoUP1008/9ZhfVbVMmTKqqrpnzx5t3bq1qqqmpqZq3bp1M8ps2rRppmeXKVNGY2JitEGDBnrkyJFM5auq9u3bV+Pi4rK0u379+vryyy+rqurcuXN1xIgRHuvYqlUrjYmJ0W7dumXUIyteffVV/fDDDzPOY2JitE2bNpqWlqb33Xefvv/++x7vzY4zZ85oWlpaprSwsDA9ePCgpqWl6fjx43X48OE53lOxYkWNj4/P9lmZPle+JOW86vcvqj4brvpSPdWti3Jd1L7jZ/QfH6zTyx5frLe8t1b/OnrGh4ZmBtigXnxXncysLgs8BjQVkVJY0UuMqo71qWf0IeqlHlGOpJ5zlu6AqVOnMm/ePAD27dvHzp07qVatGnXq1KFVq1YAtG3blr1793rMHx4enlFeZGQk4eHhbNq0iUOHDtG6dWvCw8M9Nmu+++47hg4dmqHyWKlSpYxrS5dm3fI+e/YsCQkJPPzww4C1WWVYWBjz589nyZIlnDx5krvuuotevXqRnJzM7t27GTZsGO+99x6dOnXizJkz3HvvvZQsWZLu3btnbAywYsUKPvroo4znTJgwgUmTJiEiNG7cONuIaPLkyRw7dozw8HCOHTvG3XffnRHx1K1bl3Xr1hEZGZnxfzt9+jS9e/emRIkSdOjQgbfffjtTeVndU65cOQLS33hoK8wfAwc2Q7PB0PdlKBOe831upKUp//v1b6Ys3QbAcwObcWuHCIrlczMsK5w4oiv0wlyfUKzJh9krigUYr0fN+k7J/vprzexmmRsVasOdS3Jt3/fff8/KlStZu3Ytl1xyCd27d8/YF6tUqQu7JYSEhJCYmJhtflfuvvtuZsyYwcGDBxk5cmSu7fPE1q1badu2bUYTKioqimbNmjFw4EAGDhzIiRMneOyxx+jVqxfbtm2jffv2HD9+PCP/3LlzGTJkCP379+fmm29m+PDhnD17lvj4eGrUqAHAL7/8wvLly9m0aRP33XcfSUlJNG/enOTkZJ5++mnOnj1LWloaU6dO5ZdffuGLL75g8ODBfPHFF4wePTrDCQFs2rSJqlUvzIWNjo7m6quvZvny5R7rmNU9+d4sS02Bn163+oNCK8BNn0CTAbkqat/xs/x7dhRrdx+jS/3KTBncnFoVcy8B62u8mdAodpS1Oz1NVZOA3+xXRh7/mZk7fBYRXT0x8/od8Mk2KgkJCVSsWJFLLrmE7du3s27dOp/kHzRoEBMnTiQ5OTlDSK1cuXJZRkU9e/Zk0KBBPPLII4SHh3P8+PFMUVFWREdHZ0RrYDmiAQMufEGee+65jNGlzZs307lzZ/7xj38waNAgvvvuO2JjYzO+1OnOadWqVfTo0SOjjCeffJJFixZxzTXXAGREd9OmTSMxMZGwsDD27NkDQKNGjejWrRsPPfQQR48e5f77789kb/Xq1TOdR0VF0bJly2zrmNU9+eqIDm+3oqD9m6DpILjuFWuZhhe4TjqsHhbKFXXDWR5zkGIiTL6xObe0r13g9mDzprN6lYg8ICKZpieLSEkR6SkiHwN3+Me8AkKLm6D/VCsCQqy//afmeTeDPn36kJKSQuPGjRk3blym0aK85C9ZsiQ9evTgpptuyviih4eHc+WVV9KsWbOMzmqApk2bMn78eLp160bLli155JFHMq5dd9117N+//6Ly3R1RTEwMzZo1Q1V5/PHH6du3L23atAEsR9SsWTMaNmzIiy++yE033UStWrWIjbWWIqSlWVPTli1bRp8+1uqflStXcv78+QwnBHDppZdy+vRpNm7cyJQpU3jmmWf4+OOPAfj9999p2bJlxt+ciI6OpkWLFjnmc79n2rRpREZGEhkZSadOfpIITk2BNf8H711lbeMzdIb1cuCEnpgbTVx8Igrsj09i7m9x1K50CSse7sqwDhEFzgmBF8JodjNsJNYSjzpAPBAKhABfA++o6iY/25krXJpmo3bu3JnpWoEXsMoDaWlptGnThlmzZtGgQYN8e+7UqVP5+OOPad++Pa1atWLMmKznnZ45c4b777+f0NBQunTpwvDhw2nTpg2//PILJUpkPzd20aJFfP7559SuXZuePXvSp08fXn/9dbp06cKPP/5Ily5daNcuRx0uv5Gnz9WRHdbs6LiN0PgGuP7/oGwVR0V41lcP5adxV+fOrjzgrTCaI4VGewZ1ZSBRVePzYF++UpQUGrdu3Uq/fv0YNGgQr776aqDNKXI4+ly56gaFlrd2VQ2tANe/Yg3N5yJyqTNuCVl9owXYM+V6x+XlFb8oNKo1g/pArq0y+J0mTZqwe/funDMaAou7blBSgrVGrMd4a2TMIWlpyoyf92bphMA3Uh3+JKh38fDpzGqDwZd8++zFukGaao2SOWTv0TPcMm0dkxZvpXG1coQW949Uhz8JakekRWitmaEQcXSXZ90gT+lZkJamfPTTHvq8sZptB0/y8pAWLH3oKqYMbkHNsNII1t57k29s7pcV877E66aZiLyoqo/nlGYwGDyQlmrtpPrtJKxemywaUl7qBv117AxjZ0fx657jdG9Uhck3Nqd6Bav55S+pDn/iJCK6Nos054tcChAFcOqToRCT7efp2J8w43pY8STU7Q59X7Lmornixdy0tDRlxk976PP6GrYdOMlLQ1rw0Yj2GU6osOLNhMZ7gHuBeiIS5XKpHPBT1ncVDLKbWR0aGpqxJKAgzqswFC5UlWPHjhEaGpr5Qloa/DoNVj4DISVh4LvQ8hZrRKx0mCPdoOyioMKON/OIymNJgEwGXIXtT6kljlbgyWr4Pjk5mdjY2CyXSBgMuSE0NJRatWpdmAt1fDcsuB/++gka9IL+b0D5Go7LTUtTPlm7lxeX76B4MeGp/k0Y2rZWofgB9eXw/duqepuIrFPVv3xgW4GgRIkS1KlTJ9BmGIKRtDRY/wGsfBqKFYcBb0Or4bmaF/TXsTP8e3YUv+w5TreGVZgyOHiiIFe8cURtRaQGMFJEPsHqZcugsERFBkO+cGKvFQXtXQP1r7GXBjnvOHaPgl4a0qLQREG5wRtH9C7wLVAXe5GrC2qnGwxFm7Q02Dgdvp4IUgxueBNa35arKOjvY2cZO3tz0EdBruToiFR1KjBVRN5R1XvzwSaDoXBx4i9YeD/sWQ11e1hOKKx2zve5kZamfLruL6Ys225FQYNbMLRd8EZBrngzavajqnYBbhORf7hfV1XnW0fmE2YXD4NfUYWNH8HXT1nn/V6HtiO8joJc5Tqqli9F2VLF+fPIGbo1tEbECvqyDF/iaNFrYSWrUTODIU/E/w0LH4Dd30OdbjDgLQjzfiPPdLmOxOTMWtm3tK/N5BubB00U5JdFrwZDkUcVfvsYVkwATbOkOtqNdNwX9PKKHRc5IYA1O48GjRNygpMlHqWAwUCk632qOsn3ZhkMBQRXqY5y1awtew5vhcirrCioYqTjItPSNEvNILC29ymKOImIFgAJwEYg76rxBkNBx12q49QB69VyGAx4B4o5XzOerh3tiaLUL+SKE0dUS1X7+M0Sg6Gg8e2ki6U6APb+6NgJpe+gMXnpNoqJcHO7WizcvJ/E5Ixd3AuFXIe/cOKIfhaR5qoa7TdrDIaCgmrWO7eAI6kOgNgTZ3l8ThQ/7bJ20HhxiCXT0ale5YxRsxphpRnbu1GhWzXvK7wZvo/GmrhYHLhTRHZjNc0EUFV1pkLuA0RkIHA9UB74UFW/zm8bDEHMyQOw6CHP172U6lBVPv/1b15YYu0j9sKg5gzrcGEHjcIo1+EvvImI+vnygSIy3S7zsKo2c0nvA7yBJcr/gap63GxMVecD80WkIvAKloi/wZA3VCHqK1j2b0g5Dy1uhm0Lc7WNVFx8IuPmRLFm51GurB/Oi4NbFKh9xAoa3sys/gtARIYCy1X1lIhMANoA/wGcLoSdAbwFfJKeICIhwNtYmkexwHoRWYjllCa73T9SVQ/bxxPs+wyGvHHqICz6F/yxDGpfAQPfgfB6EHWNI6kOVeXL9ft4fsk20lR5bmAzhncsmFv4FCSc9BE9paqzRKQLcA3wMtY6tI5OHqiqq0Uk0i25A7DLZSfZL4EBqjqZLCIysd7VKcAyVXVf/5aeZzQwGiAiwvuJZoYihipEz4KlYyElCXo9D1fcA8Ws/eBocZPX+9ftj0/kcTsK6lzPioJqVzJRkDc4cUTps6+uB6ap6hIRec5HdtQEXHsGY8newT2A5QwriEh9VX3XPYOqTgOmgTWz2kd2GoKJU4dgySOwfTHU6mBFQZWd7wOnqszcsI/nFm8jVZX/DGjK8I6XFYg95QsLThxRnIi8h9V8etGe4BgQ8f30hbg55TNrzQxZogoxc2DpY3D+LFz7H+h034UoyAEHEhIZNyeaH/44whV1K/HS4JZEhJsoyClOHNFNQB/gFVWNF5HqwNgc7vGWOMB1uXItOy1PqOoiYFG7du1G5bUsQ5Bw+ggseRi2LYKa7WDgf6FKQ8fFqCqzNsbyn0VbSUlTnr2hKbddYaKg3OK1I1LVs8Bcl/MD+G6zxfVAAxGpg+WAbgFuzWuhJiIyZCJmrhUFnTsF1zwLne6HEOfLLQ8mJPHE3ChW7ThChzqVeHlICy4LL+MHg4sO+b7oVUS+ALoDlUUkFnhaVT8UkfuBFVgjZdNVdUten2UiIgMAZ47Ckkdh63yo0caKgqpe7tWtrlIdNcJC6dqwCoujDpCcmsbT/ZtwR6dIEwX5gKCWAXGJiEbt3Lkz0OYYAsHWBbD4ETh3Ero/AZ0f9DoK8iTVUbdyGaaPaE9kZRMF5YS3MiBedzaLxT9EZKJ9HiEiHfJipL8xO70WYc4cg1l3wszbrfk/o3+Aqx5x1BTzJNWRlJJqnJCPcdI0ewdIA3oCk4BTwBygvR/s8gmmj6gI4SrXcUkla2Z0ShL0nABX/gtCSjgu0pMkx4F4swWVr3Ey/N5RVe8DkgBU9QRQ0i9W+QgTERUR0uU6EvYBCmePwfnT0ONJ6DrWsRNSVeZtinXbr+YCRVWqw584cUTJ9lIMBRCRKlgRksEQWLKU61DYMN1xUYdPJTH60408/NVmIiqWplTxzF+RoizV4U+cOKKpwDygqog8D/wIvOAXq3yEiPQXkWkJCQmBNsXgL84e94lch6qy4Pc4er22mh/+OML46xrz3WM9eHGwJdkhQM2w0ky+sblZMe8HHI2aicjlwNX26bequt0vVvkYI54fpOxYbsl1nD6Y9fUKteHhmByLOXLqHBPmR7NiyyFaR4Tx8pCW1K9a1sfGFk38MWo2FIhT1beBSsALItImDzYaDLkj8QTMGwNf3AyXhEPPpyx5Dle8kOtQVRZt3k+v135g1Y4jPNH3cmaP6WycUADI7er7nlg6QP/F4er7/MSMmgUhf3xtdUyfPgxd/211RhcvaW3l40Cu4+jpczw1P4ZlMQdpWTuMV4e2oH7VcvlYEYMrXjfNRGSTqrYWkclAtKp+np7mXxPzjmmaBQGJ8bBiPPz+GVRpDIP+CzVy99FbEnWApxbEcDophX9d24DRV9WleEhA1m8HPf7Y16zArL43FDF2rrQ2Mzx9CK56FLo9DsVLOS7m2OlzTFywhSXRB2hRqwKvDG1Jw0tNFFQQKCir7w2Gi0lKsKKgTZ9Clcvhls+gZttcFbUs+gAT5sdwMimZsb0b8c+uJgoqSBSU1fd+wfQRFWJ2fQsLH4RT+6HLw9BtHJQIdVzM8TPneXrhFhZt3k/zmhX4fOgVNKpmoqCChqPV97ZYfQMg4xOhqqt9bZSvMKvvCyFJJ+HrCda2zpUbwl3fQK0cuxiyZHnMQSbMjyYhMZlHr23ImO71KGGioAKJky2n7wYewhIt+x24AliLNYJmMOSdP1dZfUEn46xV8j3Gex0Fucp1VKsQSvUKofz2dzxNa5Tn07s60rh6eT8bb8gLTiKih7AWuK5T1R725MYCPbPaUEg4dwq+mWgtyQivDyNXQG3vhR3c5ToOJCRxICGJPk2r8eatrU0UVAhw4oiSVDVJRBCRUqq6XUTMohtD3tizGhbcB/H7LMUKImirAAAaOUlEQVTEnhMunpyYA57kOqLjEowTKiQ4cUSxIhIGzAe+EZETON/TzFCUcZXqKF8DwhvAnu+hUj0YuRwirshVsXEe5Do8yXgYCh5ORs0G2YfPiMgqoAKwzC9W+QgzalaASJfqSF8lfzLOetW7Fm7+BEo63/kiITGZSYu2erxu5DoKD07WmpUSkVtF5EmgG9AKeMJvlvkAo0dUgMhSqgM4uj1XTmjV9sP0eu0H5v8eR+8mVQktYeQ6CjNOmmYLgARgI3DOP+YYghZPkhwOpDrAioKeW7yVWRtjaXhpWT64vT3Na1VwE7kvzdjejYxcRyHCiSOqpap9/GaJITg5f9aKhvCwprFCLa+L+n7HYcbNiebI6XPc16MeD17dgFLFrU0RB7auaRxPIcaJI/pZRJqrarTfrDEEF3+thQX3wvHdUKcH7FsHKS7NMy+kOgBOJiXz/OJtfLVhHw2qluW929rSsnaYHw035Dc5OiIRicb6OSsO3Ckiu7GaZgKoqrbwr4mGQkdyInz3HKx9G8Jqwx2LoE7XzKNmXkh1AKz+4wjj5kRx8GQS93Svx0NXNyC0hPOtoQ0FG28ion5+t8IQPOz7FebfA8d2Qbu74NpJUMoWGmtxU46OJ51TScm8sHQbX/y6j3pVyjDnns60jqjoR8MNgSRHR6SqfwGISChwL9AFK0L6EUsYLV8RkcZYs7wrY8nV5rsNhixIToRVz1tRUPmacPsCqNs9V0X9uPMoj8+J4kBCIv/sWpeHr21ooqAgx0kf0SdYe5m9aZ/fCnwKDPW2ABGZjhVhHVbVZi7pfYA3sLab/kBVp3gqQ1W3AWNEpJhtk3FEgWbfeqsv6Ogf0HYEXPsfCHW+tuv0uRQmL93G/375m7qVyzBrTGfaXmaioKKAE0fUTFWbuJyvEhHPs8myZgbwFpYDAcDeouhtLMG1WGC9iCzEckqT3e4fqaqHReQG4B4sR2gIFMlJ8P0L8PObUK4G3DYP6uVuDfTPfx7l37OjiItPZNRVdXi0VyMTBRUhnDii30TkClVdByAiHQFH+ququlpEIt2SOwC7VHW3Xe6XwABVnYyH/ilVXQgsFJElwOdZ5RGR0cBogIiICCdmGrwhbiPMuweO7oA2t0Ov53MVBZ05l8KUZdv5dN1f1Klchln/7ES7yEp+MNhQkHHiiNpiDeH/bZ9HADvSR9XyMHpWE3DdmCqWbAT5RaQ7cCNQCljqKZ+qTgOmgaVZnUvbDO6knIPvp8BPr0O56vCPOVD/Gq9udZ90OKh1DRZs3k/siUTu6lKHx3o1onRJEwUVRZw4ogIxmVFVvwe+9yavWWvmY+J+g/n3wpFt0Pof0PsFCPVu+Yy7VEdcfCJvrfqT8DIl+Gp0JzrUMVFQUcbJold/rbSPA2q7nNey0/KMUWj0ESnn4IeX4MfXoGxVuHUWNOzlqAhPUh0li4cYJ2RwJhXrJ9YDDUSkDpYDugVrRC7PmIgol7hOPCxbFSTE0o5ueSv0eQFKOx/J8iTJcTAhKa/WGoKAfFWNEpEvsORlG4lIrIjcpaopwP3ACmAbMFNVt/jieWb1fS5Il+tI2AeotYXPqf3Q6QFrL7FcOKENe49TrJhkec1IdRjAy4hIRARr0eu+HDNng6oO85C+lGw6nnOLiYhygSe5jq3zofdzjopKSk7llRU7+PCnPYSVLsGZ86mcT0nLuG6kOgzpeBURqbUdrM8dhb8xEZFDUpPtSCgLHMp1bPzrBNe9sYYPftzD8I4R/Ph4T14a3IKaYaURoGZYaSbf2NysmDcAzucRtVfV9X6zxseYiMgBB2OsNWKe8FKuIyk5lf/75g8+WLOb6hVK87+7O3Jl/cqAkeoweMZJH1FHYK2I/CkiUSISLSJR/jLMF5iIyAtSk60RsWnd4dQBuOLei8XrvZTr2PT3Ca6fuoZpq3dzc/sIlv/rqgwnZDBkh5OIqLffrDAEhkNbYf4YOLAZmg2Gvi9DmXCo0dqRXEdSciqvrfyD91fvplr5UD4Z2YGuDavkY0UMhR1H84iy2umVAryTh2maeSA1xZoZ/f0Ua0LiTZ9AkwEXrjuQ69i8L55HZ21m1+HT3NK+Nk9e35jyoSX8ZLghWAnqnV7NhMYsOLzdioL2b4Kmg+C6V6CM8+bTuZRU3li5k3d/+JNLy4fy8cgOdDNRkCGXmJ1eiwqpKbD2TVj1ApQqB0NnWI4oF0TFxvPYrM38ceg0N7WrxYR+TUwUZMgTQb3Tq2ma2RzZYY2IxW2ExjfA9f8HZZ1HL+dSUnnz213894c/qVy2JB+NaE+Py6v6wWBDUSOod3ot8k2ztFRY+xZ89zyULANDpkPTG0GynuWcHTFxCTw6czM7Dp1icJtaTOzXhAqXmCjI4BvyutPrcr9YZcg7R3daUVDseri8H/SzF6x6gatcR/UKoTSvWYGV2w8TXqYkH97RjqsbX+pn4w1FDSed1QIMB+qq6iQRicDa7fVXfxlnyAVpqbDuHWsXjeKhcOMH0HyI11GQu1zH/oQk9ick0e6yMD68o4OJggx+wUnT7B0gDWuUbBKWfvUcrA7sAkmR6yM6usvSjt73CzS63oqCyjmLXjzJdRxIOGeckMFvOHFEHVW1jYhsAlDVEyJS0k92+YSg7iPKtEdYTYjoDNsWWlHQoGnWPKBc9AV5kuvwlG4w+AInjijZFrpXABGpghUhGfKbdKmO9FXyCbEQPROqtYBbZ0L56o6LTE5N47/f/+lpY2gj12HwK04c0VRgHnCpiDwPDAEm+MUqQ/Z4kupIPJErJ7Tj4CkenfU7MXEnaV07jG0HT5KUbOQ6DPmHk1Gz/4nIRuBqO2mgvceYIb/xJMnhUKojJTWN91bv5vWVf1A+tAT/Hd6Gvs2rXyRyP7Z3I7Nq3uBXnIyalQLaYA3bFweGigiqOslfxhncSEuDDR96vu6lVAfAzkOneHTWZqJiE7i+RXUm3dCU8LKlACPXYch/nDTNFgAJwEbgnH/M8S1BNWp2Yi8suB/2roGqTeH4n5DiovfspVRHSmoa76/Zw2vf/EHZ0OK8fWsbrm/hvDlnMPgSJ46olqoWiC2FvCUoRs3S0mDjdPh6IkgxuOFNaH0bRM9yJNUBsOvwKR6dFcXmffH0bVaN/wxsRmU7CjIYAokTR/SziDRX1Wi/WWPITPzfVhS05weo28NyQmH2zksOpDpS05QP1uzm1W/+oEzJEN4c1pp+LaojuRjeNxj8gRNH1AUYISJ7sJpmQt52eDV4QhU2fgRfP2Wd93sd2o7I1bygP4+c5rFZm9n0dzy9m17KcwObU6WciYIMBQsnjqiv36wwXCB+Hyx8AHavgjrdYMBbEBbhuJjUNGX6j3t45esdhJYI4Y1bWnFDyxomCjIUSArCTq8GsKKg3z6BFeNB0yypjnYjcxUF7T5ymrGzo9j41wmuaXwpLwxqRtXyoTnfaDAEiBwdkYj8qKpdROQUZJp4m940K+8364oKCXFWFPTntxB5lRUFVYx0XExamvLRz3t5afl2QkuE8NrNLRnYqqaJggwFnhwdkap2sf+W87853iEiZYAfgGdUdXGg7ck1qrDpM1jxJKSlWLKt7e6CYjlvruI+6XBE50i+2XqIX/ce5+rLq/LCjc251ERBhkKCkz6iPCMi04F+wGFVbeaS3gd4AwgBPlDVKTkU9Tgw02+G5gcn98PCB2HXN3BZFysKqlTHq1vdpTri4hN5fuk2ShUXXh3akhvbmCjIULjwpmmW3iTL6pPttGk2A3gL+MSl/BDgbeBaIBZYLyILsZzSZLf7RwItga1k3kmk8KAKm7+AZeMg9Tz0fQnaj/IqCkrHk1RHxUtKMbit97OrDYaCgjdNM581yVR1tYhEuiV3AHap6m4AEfkSGKCqk7Gip0yISHegDNAESBSRpap6kQqAiIwGRgNERDgfdfIZrnId5arDJeFwKBoiOsGAtyG8nuMiPUlyHDqZlGW6wVDQcbrWbDAQ6XqfD9aa1QRcN1yPxdpVNktUdbxtzwjgaFZOyM43DZgG0K5dO0/qFv7FXa7j1H7r1fwmGPSeoygonX3Hz1KieDHOp1xcbSPVYSisFNq1Zqo6I6c8AV9r5kmu4++1jp1QWpryv1//ZvLSbYgqJUKE5NQL/tVIdRgKMwVhrVkcUNv1OXZangnoWjNVSNiX9TWHch2xJ87y+Jwoftp1jKsaVGbK4Bas33PcSHUYgoaCsNZsPdBAROpgOaBbgFt9UXDAIqLTh2Hxw56veynXoap8/uvfvLDEkn16YVBzhnWojYhQ00h1GIIIJ+2DLsBGEdkhIlEiEi0iUU4eJiJfYG1T3UhEYkXkLlVNAe4HVgDbgJmqusVJuZ5Q1UWqOrpChQq+KM6bB0L0bHi7I+z8BpoNseQ5XPFSriMuPpHbp//K+HkxtIoIY8XDXbm1Y4QZljcEJfm61kxVh3lIXwoszWv57uRrRHT6CCx5xBKwr9kWBv4XqjRyE7nPWa5DVflq/T6eW7KNNFWeG9iM4cYBGYIcUfV+QElEWgJX2adrVHWzX6zyMe3atdMNGzb47wFb5sGSR+HcKejxJHR6AEKczxU9kJDI43OiWf3HETrVDeelIS2oXekSPxhsMOQPIrJRVdvllM/J8P1DwChgrp30mYhMU9U3c2mj3/F7RHTmqOWAts6HGq2tKKhqY8fFqCqzNsbyn0VbSUlT/jOgKcM7XkaxYiYKMhQNvI6I7P6gTqp6xj4vA6wtDHpEfomIti6AxY9AUgJ0HwdX/itXUdDBhCSemBvFqh1H6FinEi8PaUlEuImCDMGBzyMirCUerusKUsl62Udwc+YYLH0MtsyF6q3gjkVwaRPHxagqc36L49lFW0hOTeOZ/k24vVOkiYIMRRInjugj4BcRmYflgAYA0/1ilY/wedNs2yJrWD4xHnpMgC7/ghDn2zAfOpnEE3Oj+W77YdpHVuTlIS2JrFzGNzYaDIUQp53VbbCG8RX4UVU3+cswX5LnptnZ47Ds35ZgfbXmMPBdqNYs5/twl+sIpXujKizafIDzqWmM7X05d3Y2UZAhePFHZ3Up4HKgrH1ffxHpH/T7mm1fCov/BWePQfcn4apHvI6CLpbrSOJ/v+wjMvwSpo9oT90qZf1pucFQaCi0a828IU9Ns7PHYfk4iPoKLm0Gw2dDdWf98p7kOs6nphknZDC4UBDWmvkNr9eauU86bDLAmiF99ih0exyuegyKl3T8fE9yHQfijVyHweBKQVhrFljcpToS9sHat6BcTbj7W6jRynGRqsqiqAOIWKs+3DFyHQZDZsy+Zp6kOopJrpzQ0dPneGp+DMtiDlK7UmkOnzzHORftICPXYTBcTFDva+ZVH5EnSY4E50okS6IO8NSCGE4npfB4n8sZdVUdFkcdMHIdBkMOOBq+L6xkO3z/WrOsdYMq1IaHY7wq/9jpc0xcsIUl0QdoUasCrwxtScNLC8ymJwZDwPDHzOrg5OqJmfuIwGupDoBl0QeYMD+Gk0nJjO3diH92rUvxEOcSsAZDUcY4onRJDgdSHQDHz5xn4oIYFkcdoHnNCnw+9AoaVTNRkMGQG4wjAsvp5OB4XFkec5AJ86NJSEzm0WsbMqZ7PUqYKMhgyDXGETngxJnzPLNoCwt+30+T6uX59K6ONK5udtw2GPJKUDsiXy56/WbrIZ6cF82JM+d5+JqG3NvDREEGg68Iakfki1084s+e59lFW5m3KY7Lq5Vjxp3taVojnzSwDYYiQlA7orzy7bZDPDE3muNnzvPg1Q24v0d9ShY3UZDB4GuMI8qChMRkJi3aypzfYrm8Wjmmj2hPs5omCjIY/IVxRGTWDKpUpiTJqamcOZ/GAz3r80DPBiYKMhj8TJF3RO6aQcfOnEeAh69twINXNwyscQZDEaHQ/dSLSHcRWSMi74pI97yWl5VmkAJfrXe2LbTBYMg9+eqIRGS6iBwWkRi39D72DrK7RGRcDsUocBoIBfLsLTxpBnlKNxgMvie/m2YzgLeAT9ITRCQEeBu4FsuxrBeRhUAIMNnt/pFYGzv+ICKXAv8HDM+LQTXCShOXhdMxmkEGQ/6RrxGRqq4GjrsldwB2qepuVT0PfAkMUNVoVe3n9jqsquniPieAUnm1aWzvRpQuEZIpzWgGGQz5S0HorK4JuOpwxAIdPWUWkRuB3kAYVnTlKd9oYDRARESEx4enawMZzSCDIXAUBEfkCFWdy4Vtr7PLNw2YBpYeUXZ5B7auaRyPwRBACsKoWRxQ2+W8lp2WZ0Skv4hMS0hI8EVxBoPBTxQER7QeaCAidUSkJHALsNAXBavqIlUdXaGCmRVtMBRk8rVpJiJfAN2ByiISCzytqh+KyP3ACqyRsumqusVHz+sP9AdOishhrH3Z3KmQRbp7mut5TseVgaN5MDsre7y9XpDqklM9csqTnd05nacfu6YFqi5O3xP3c/e6+PvzlV2e3Hy+LvPKKlUtEi9gmrfp7mmu5zkdAxv8Yac31wtSXXKqh9O6ODl3sd81LSB1cfqe5FQXf3++fFkXb56V/ioITbP8YpGDdPe0RQ6P80JO5WR3vSDVxZsynNTFyfkiD3lyS17q4vQ9cT8vzHXx2uYisYtHfiIiG9SLXQsKA6YuBY9gqYc7RSkiyi+mBdoAH2LqUvAIlnpkwkREBoMh4JiIyGAwBBzjiAwGQ8AxjshgMAQc44jyERGpKyIfisjsQNuSG0SkjIh8LCLvi0ie5FcCSWF/H1wRkYH2+/GViPQKtD25xTgiL/GFqJtaUid3+ddSZzis143AbFUdBdyQ78Zmg5N6FMT3wRWHdZlvvx9jgJsDYa8vMI7Ie2YAfVwTXETd+gJNgGEi0kREmovIYrdX1fw32Stm4GW9sBYkp0u2ZNbXDTwz8L4eBZ0ZOK/LBPt6oaTQyYAEClVdLSKRbskZom4AIpIu6jYZ6Je/FuYOJ/XC0oqqBfxOAfsRc1iPrflrnTOc1EVEtgFTgGWq+lu+GupDCtSHqRCSlaibR2EjEQkXkXeB1iLyhL+NywOe6jUXGCwi/8V3Sw78SZb1KETvgyue3pMHgGuAISIyJhCG+QITEeUjqnoMqy1fKFHVM8CdgbYjrxT298EVVZ0KTA20HXnFRER5w2+ibgEmWOoVLPWA4KrLRRhHlDf8JuoWYIKlXsFSDwiuulyEcUReYou6rQUaiUisiNylqilAuqjbNmCm+kjULb8IlnoFSz0guOriLWbRq8FgCDgmIjIYDAHHOCKDwRBwjCMyGAwBxzgig8EQcIwjMhgMAcc4IoPBEHCMIzLkCyLys4f0GSIyJL/t8QUiMtvWNnpIRF53SX9PRFa6nD8gIlNFpKSIrBYRs7TKDeOIDPmCqnYOtA25RSyKuaU1BULs1fA/Aa71awlUsKU7sK/9rKrngW8pxLpB/sI4oiBBRMaKyIP28Wsi8p193FNE/mcf/1dENojIFhF51k7rIyKzXMrpLiKL7eNeIrJWRH4TkVkiUtZOv05EtovIRvuXPj3/MyLymEtZMelyFiJy2v4rIvKWLfC1Eqjqkr+tiPxgl7tCRKpnUc8ZIvKuXY8/RKSfnR4qIh+JSLSIbBKRHnb6EhFpYR9vEpGJ9vEkERnl8r9bLyJRLv+XSNvGT4AYMq/zAhgOLLCPfwcaikhpEakAJNppze3rnbGcFcB8+16DC8YRBQ9rgKvs43ZAWREpYaetttPH25vztQC62V/QlUBHESlj57kZ+FJEKmOJbV2jqm2ADcAjIhIKvAf0VdW2QBWHdg4CGmGJe92OHUnYtr4JDLHLnQ4876GMSCx9nuuBd22b7gNUVZsDw4CP7fQ1wFW2g0gBrrTLuApYLZa8agO7vFZAWxHpaudpALyjqk1V9S83G64ENmI9NAXYBLQHrgB+AdYBnUWkJtYKhnQJjxg7n8EF44iCh41YX6LywDmstUrtsL5wa+w8N4nIb1hfmqZAE/tLtBzob/ddXI/1S38FlrP4SUR+B+4ALgMuB3ar6h67zC8c2tkV+EJVU1V1P/Cdnd4IaAZ8Yz9vAtYK86yYqappqroT2G3b1AX4DEBVtwN/AQ3tunfFchxLsBz0JUAdVd0B9LJfm4Df7LIa2M/5S1XXebChOnDE5fxnLKfaGet/v9blPKN/TFVTgfMiUi67f1JRw3SaBQmqmiwie4ARWB/8KKAHUB/YJiJ1gMeA9qp6QkRmAKH27V9iLag8DmxQ1VMiIsA3qjrM9Tki0iobM1LI/OMW6iljFgiwRVU7eZHXfYFkdgsm12M55N3AN0BlYBR2NGM/d7KqvpfJGKtJeSabchPJXL+fsDSOQrEkW49gOfIjuDgim1JAUjZlFzlMRBRcrMFyNqvt4zHAJrVWNpfH+mIliMilWNrH6fwAtMH6gn5pp60DrhSR+pCxg0dDYAdQVy5Imbp2vO61y0FE2gB1srBxNXCziITYfUA97PQdQBUR6WTfX8LuEM6KoSJSTETqAXXte9dg973YdkYAO+wO4n3AUKwoxfV/BNZq9pEu/V81xTt98W1YTj6dtVhRZBVVPWz/z49gSdOm9w8hIuHAUVVN9uIZRQbjiIKLNVhNhrWqegjrV3cNgKpuxmp+bAc+x+XLYTcXFmM5p8V22hGs6OoLEYnC+qJdrqqJwL3AchHZCJwCEuyi5gCVRGQLVoT1RxY2zgN2YulGf2KXi+0whgAvishmrM5eTyNtfwO/AsuAMaqaBLwDFBORaOArYISqnnP5vxy2bV+D1eRL/798bf8/1tr3zga8aTYtAbqnn6jqCSzH4yrNsRarM36zS1oP+16DC0YGxOAYESmrqqft5tvbwE5VfS2fnj0DWKyqAd2TTERKA6uAK21H7u19c4FxqpqVky6ymIjIkBtG2R3KW4AKWKNoRQo7unqabDZLcEcsZcX5xgldjImIDAZDwDERkcFgCDjGERkMhoBjHJHBYAg4xhEZDIaAYxyRwWAIOMYRGQyGgPP/Z+DWxrCKQl0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# compute the analytical form of the analytical phase shift from SPM\n",
    "\n",
    "# get effective area\n",
    "field_profile = simulation.fields['Ez'][Nx//2,:]\n",
    "width_eff = np.square(np.sum(np.square(np.abs(field_profile)))) / np.sum(np.power(np.abs(field_profile), 4))\n",
    "width = dl*width_eff\n",
    "height = width*100\n",
    "Aeff = width*height # Assume square wg if extrapolated to 3D\n",
    "\n",
    "# get the analytical self phase\n",
    "n2 = 3*chi3/simulation.L0/simulation.L0/(C_0/simulation.L0)/np.sqrt(eps_m)/(EPSILON_0*simulation.L0)\n",
    "L = dl*(Nx-200)\n",
    "gamma_spm = (omega/3e8*simulation.L0)*n2/Aeff\n",
    "\n",
    "# plot the nonlinear phase shift vs power\n",
    "plt.figure(figsize=(4,3))\n",
    "plt.loglog(pwr_vec*height, -np.unwrap(np.angle(T_vec)-np.angle(T_linear))/np.pi, \"-o\", label=\"fdfd\")\n",
    "plt.loglog(pwr_vec*height, (pwr_vec*height)*L*gamma_spm/np.pi, \"-o\", label=r\"analytic: $n_2k_0/A_{eff}\\cdot P \\cdot L$\")\n",
    "plt.xlabel(\"waveguide power (W)\")\n",
    "plt.ylabel(\"nonlinear phase shift ($\\pi$)\")\n",
    "plt.title(\"Waveguide Kerr SPM\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
